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INTRODUCTION: 

The  following  statement  from  our  1995  IDEA  grant  proposal  has  held  up  over  the 
intervening  4  years: 

“Early  detection  of  breast  cancer  is  critical  to  successful  treatment.  Unfortunately, 
conventional  x-ray  mammography  has  severe  limitations  in  early  tumor  detection,  in  particular  for 
women  with  radiographically  dense  breasts.  For  this  reason,  nuclear  medicine  techniques  are 
beginning  to  be  applied  to  diagnostic  breast  imaging.  Among  nuclear  medicine  techniques, 
Positron  Emission  Tomography  (PET)  has  the  greastest  potential  for  high  sensitivity  to  small 
tumors,  and  good  uptake  ratios  have  been  demonstrated  for  breast  tumors  with  the  standard  PET 
radiotracer  Fluorodeoxyglucose  (FDG).  The  key  PET  detector  attributes  for  this  application  are 
high  acceptance  and  good  efficiency  at  high  rates,  fine  spatial  resolution,  and  good  scatter 
rejection.” 

For  these  reasons,  we  proposed  to  construct  and  test  modular  components  of  a  novel 
prototype  PET  detector  module  which  was  to  produce  very  high-resolution  images  in  all  three 
dimensions,  with  large  acceptance  and  high  rate  capability.  We  did  build  and  test  such  modules, 
and  their  performance  was  in  general  in  reasonably  good  agreement  with  expectations.  Our 
original  plan  was  to  design  modules  as  parts  of  a  unit  to  be  mounted  on  a  standard  mammography 
stand.  The  prototype  devices  as  constructed  were  however  not  at  all  ergonomically  convenient 
when  deployed  on  a  mammography  stand  (although  we  did  acquire  and  mechanically  modified  a 
mammography  gantry  for  this  purpose),  despite  our  best  efforts  within  our  budget.  Subsequently, 
we  have  continued  to  work  on  positron  emission  mammography  through  research  support  from 
other  sources  after  our  US  Army  grant  effort  finished,  and  recently  we  have  developed  a  system 
which  is  much  more  compact  and  is  compatible  with  x-ray  systems.  In  what  follows  we  will 
discuss  in  detail  our  experience  in  the  US  Army-funded  development  program  which  established 
the  feasibility  of  the  detector  design  which  are  continuing  to  carry  forward,  as  well  as  the  important 
data  acquisition  electronics  and  image  reconstruction  techniques  which  we  developed  with  (in  part) 
US  Army  support.  The  categories  below  will  be  reproduced  in  the  Table  of  Contents  and  in  the 
Body  of  this  Final  Report. 

Detector  Construction: 


We  proposed  to  construct  our  prototype  detector  modules  using  newly-developed  LSO 
(Lutetium  OxyorthoSilicate)  scintillator  crystals  read  out  through  wavelength-shifting  fiber  optics. 
As  it  happened,  the  availability  of  LSO  scintillator  was  very  much  more  limited  than  was  projected 
by  its  sole  source  (CTI,  Inc.)  at  the  time  of  our  proposal.  We  were  eventually  able  to  procure 
some  100  cc’s  of  this  material,  and  were  able  to  characterize  devices  built  from  it  with  some 
accuracy.  This  was  much  less  material  that  would  have  been  needed  to  make  efficient  large-area 
detectors,  however.  The  sole  supplier  of  LSO  subsequently  restricted  its  distribution  to  “internal 
use  only”  within  there  organization  for  the  next  two  years.  Consequently,  much  of  our  subsequent 
efforts  were  carried  out  using  CsI(Na)  scintillator  rather  than  LSO.  While  CsI(Na)  is  readily 
available  and  a  factor  of  10  less  expensive  than  LSO,  it  is  also  near  10  times  slower  and  has  a 
smaller  photofraction.  Nonetheless,  it  proved  quite  a  reasonable  material  for  our  development 
effort  —  somewhat  akin  to  developing  a  sculpture  in  wood  or  plaster  before  attempting  it  in  marble. 
Details  of  our  experience  constructing  these  new  devices,  both  at  the  component  and  system  level, 
are  presented  as  a  narrative  in  the  following  section.  These  results  are  given  at  a  more  detailed  level 
than  in  the  publications  which  we  produced  and  have  enclosed  as  discussions  of  other  aspects  of 
our  US  Army-funded  research  effort.  We  have  a  included  a  quite  detailed  description  in  this  report 
both  as  an  indication  of  the  level  of  effort  required  and  as  a  guide  to  others  who  might  seek  to 
follow  in  our  footsteps  —  in  particular  to  alert  them  to  the  pitfalls  into  which  we  on  occasion  fell. 


Optoelectronics  and  Data  Acquisition  Electronics: 


Details  of  our  work  characterizing  and  selecting  an  appropriate  multianode  photomultiplier 
system  for  wavelength-shifting  fiber  readout  are  given  in  a  Conference  Proceeding  which  follows 
(“Characterization  of  a  new  multianode  PMT  for  low-level  optical  fiber  readout”).  One  aspect  of 
our  system  which  required  considerably  more  time  and  effort  than  originally  estimated  was  the 
development  of  high-rate  portable  data  acquisition  electronics.  The  final  system  which  we 
developed  is  discussed  in  a  recent  submission  to  IEEE  Transactions  in  Nuclear  Science  (“A  500K 
Event/Sec  12-Bit  ADC  System  with  High-Speed  Buffered  PCI  Interface”). 

Detector  Performance: 


Details  of  our  detector  system  performance  are  given  in  a  Conference  Proceeding 
(“Development  of  a  High-Resolution  PET  Detector  using  LSO  and  Wavelength-Shifting  Fibers” 
and  a  publication  which  has  been  accepted  in  a  refereed  journal  (“First  Results  with  High- 
Resolution  PET  Detector  Modules  using  Wavelength-Shifting  Fibers”  in  IEEE  Transactions  in 
Medical  Imaging.)  The  intrinsic  spatial  resolution  was  not  as  precise  as  had  been  originally  hoped 
(l-2mm  FWHM),  principally  because  the  contribution  to  resolution  blurring  from  Compton 
scattering  within  the  detectors  was  underestimated  in  our  original  estimates  of  anticipated  device 
performance.  The  best  we  were  able  to  achieve  in  practice  was  2-3mm  FWHM.  Furthermore,  at 
this  level  of  performance  intrinsic  spatial  resolution  is  significantly  less  important  than  high 
efficiency,  large  acceptance,  and  high  rate  capability.  Because  of  the  limited  activity  which  can  be 
delivered  to  the  patient,  and  because  of  the  limited  sensitivity  of  even  with  a  very  efficient  detector 
with  relatively  high  protofraction  (fraction  of  most  usable  events)  in  close  proximity  to  the  breast 
(very  large  acceptance),  accurate  detection  and  quantitation  of  very  small  lesions  is  limited  by  event 
statistics  rather  than  intrinsic  resolution  for  a  device  with  <3mm  FWHM  resolution. 

Image  Reconstruction: 

In  the  process  of  developing  image  reconstruction  techniques  for  use  with  our  novel  breast 
imaging  detector,  we  serendipitously  stumbled  upon  a  very  promising  3D  iterative  image 
reconstruction  method  using  Monte  Carlo  methods.  This  “Inverse  Monte  Carlo”  or  “Monte  Carlo 
System  Response  sampling”  techniques  has  turned  out  to  a  representative  of  an  entire  class  of 
procedures  for  implementing  “standard”  3D  iterative  reconstruction  algorithms  without  requiring 
explicit  calculation  of  the  system  response  matrix.  This  system  response  matrix  is  so  large  that 
even  when  allowing  for  symmetries  and  zero-suppressing  to  allow  for  its  sparseness,  it  is 
impractical  to  store  it  explicitly  even  on  disk,  much  less  in  computer  memory.  Our  “Inverse  Monte 
Carlo”  reconstruction  technique  is  therefore  of  quite  general  importance,  and  particularly  of  interest 
for  breast  imaging  because  of  the  required  detector  acquisition  geometry  in  this  case.  One 
Conference  Proceeding  (“An  Inverse  Monte  Carlo  Algorithm  for  3-D  PET  Reconstruction”)  and 
extended  abstracts  which  have  submitted  to  a  future  conference  (“Quantitative  Evaluation  of  an 
Inverse  Monte  Carlo  procedure  for  3-D  PET  Image  Reconstruction”  and  “3-D  Iterative 
Reconstruction  Procedures  for  PET  Detectors  with  Fixed  Dual-Plane  Geometry”). 

Future  Directions: 


As  noted  earlier,  we  have  obtained  support  for  follow-on  efforts  and  remain  committed  to 
developing  instrumentation  for  Positron  Emission  Mammography  (PEM) ,  both  using  the  detector 
design  whose  development  was  supported  by  this  grant,  and  using  other  detector  designs  which 
are  more  compatible  with  x-ray  mammography  gantries.  The  fiber-based  design  which  this  effort 
helped  to  develop  is  more  cost-effective  for  large-area  application  and  gives  important  depth-of- 
interaction  effort,  while  our  newer  more  compact  design  is  more  x-ray  compatible.  We  may  in  fact 
meld  the  two  designs  at  some  future  juncture.  The  3-D  iterative  reconstruction  techniques  which 
this  effort  helped  us  to  develop  have  wide  application  for  PEM  and  in  fact  to  3-D  PET  generally. 


DETECTOR  CONSTRUCTION 


The  details  of  our  detector  design  have  evolved  through  several  iterations  of  a 
construction/test  cycle,  and  have  been  informed  by  a  variety  of  component-level  tests  and  computer 
simulations.  Principle  component-level  test  and  development  results  are  indicated  below. 

1)  Crystal  geometry  and  surface  preparation: 

Our  crystal  geometry  of  112mm  x  112mm  x  2.5mm  was  based  upon  readout  through  8 
ribbons  each  on  the  x-  and  y-  surfaces  of  each  crystal,  where  each  ribbon  is  14mm  across.  This 
14mm  number  arose  from  measurements  on  early  production  units  of  the  Hamamatsu  R5900-L16 
position-sensitive  photomultipliers  which  we  use  to  read  out  the  fibers.  The  2.5mm  thickness  was 
based  on  a  target  resolution  of  2.0mm  FWHM  (Full  Width  at  Half  Maximum),  which  we  expect  to 
be  a  good  match  to  the  statistical  limits  in  image  quality  available  for  a  high-efficiency  PET 
mammography  unit  for  anticipated  clinical  doses  and  scan  times.  The  quality  of  the  scintillator 
crystal  surface  finish  is  important  in  our  application  for  two  reasons:  first,  because  we  require  a 
clear  division  between  light  which  is  totally  internally  reflected  within  the  crystal  and  that  light 
which  emerges  from  the  crystal  into  the  fiber  ribbons  (for  good  spatial  resolution),  and  second 
because  we  need  the  overall  crystal/fiber  laminated  stack  to  be  transparent  to  wavelength-shifted 
light  (for  good  energy  resolution).  We  require  such  transparency  because  in  our  device  we  have 
separated  the  functions  of  interaction  position  measurement  (through  the  fiber  readout)  from 
gamma  ray  energy  measurement  and  triggering/timing  (through  an  Anger  array  of  photomultipliers 
which  view  each  stack).  Because  both  the  scintillator  crystals  and  the  fibers  are  transparent  to  the 
wavelength-shifted  light,  photons  traversing  the  wide  range  of  optical  indices  within  the  stack  are 
scattered  rather  than  absorbed;  small  scratches  on  the  crystal  surface  from  imperfect  polishing, 
however,  absorb  light  and  can  cause  excessive  attenuation  of  wavelength-shifted  light  within  the 
laminate.  In  our  devices,  we  surround  the  laminated  stack  by  reflective  material  on  all  sides  except 
that  side  with  the  Anger  array,  and  we  are  able  to  obtain  good  energy  and  timing  resolution  as  well 
as  good  uniformity  of  response.  Because  CsI(Na)  does  not  readily  polish  to  an  optical  finish  and 
because  it  is  somewhat  hygroscopic,  considerable  effort  was  needed  in  order  to  obtain  an  adequate 
surface  finish  (using  silicon  oil  and  progressively  smaller  sizes  of  Aluminum  Oxide  powder 
abrasives).  Our  final  procedure  yielded  less  than  10-20%  variation  in  light  collection  from  one  side 
of  the  stack  versus  the  other,  which  is  readily  corrected  since  the  depth  of  interaction  in  the  stack  is 
measured. 

The  effect  of  imperfect  surface  preparationof  CsI(Na)  on  spatial  resolution  was  difficult  to 
analyze,  and  measurement  of  the  amount  of  scintillation  light  scattering  from  crystal  surface 
imperfections  (rather  than  totally  internally  reflecting)  required  a  great  deal  of  experimentation. 
Finally,  we  were  able  to  determine  that  our  imaging  resolution  for  2.5mm  thick  plates  of  CsI(Na) 
is  limited  by  Compton  scattering  and  secondary  photon  capture  within  a  single  plate.  The  true 
photofraction  (ratio  of  photocapture  events  to  all  gamma-ray  interactions)  for  5 1 1  keV  gamma  rays 
in  CsI(Na)  is  about  15%,  while  the  apparent  photofraction  (fraction  of  events  within  the  photopeak 
for  5 1 1  keV  gammas  normally  incident  on  a  2.5mm  thick  plate)  is  more  than  25%.  The  difference 
arises  from  gammas  which  Compton  scatter  at  90  degrees,  followed  by  photocapture  of  the 
secondary  gamma  rays.  These  secondary  gammas  have  energies  near  250  keV  and  a  range  of 
about  7  mm  in  CsI(Na);  computer  simulations  show  a  limiting  resolution  due  to  this  effect  close  to 
our  observed  2.5mm  FWHM,  for  two  detectors  in  coincidence.  For  a  CsI(Na)  device,  this  effect 
can  be  somewhat  lessened  with  thinner  crystals,  which  increases  the  fraction  of  “true” 
photocaptures  within  “apparent”  photocaptures  (by  decreasing  the  angular  acceptance  for 
secondaries  to  convert  in  the  same  crystal)  but  this  either  increases  the  number  of  detector  layers 
required  or  decreases  the  overall  efficiency. 


2)  Optical  coupling  of  crystals  to  fibers: 


Our  second  round  of  prototype  modules  (the  first  had  attenuation  difficulties  due  to 
imperfect  crystal  surface  treatment)  yielded  a  significantly  lower  fiber  light  yield  relative  to  the  first, 
for  reasons  which  were  not  understood  at  first.  After  some  experimentation,  we  discovered  that 
the  effect  was  due  to  a  relatively  high-index  optical  epoxy  which  we  used  to  coat  the  hygroscopic 
crystals;  once  we  determined  that  die  crystals  did  not  degrade  with  ambient  humidity  over  a  time 
scale  of  1-2  weeks  (the  time  between  cutting/polishing  and  sealing  them  into  laminated  modules) 
we  omitted  the  coating  step.  Although  the  crystals  did  not  degrade,  the  greater  optical  index 
mismatch  without  the  coating  trapped  more  light  within  the  crystal  and  decreased  our  fiber  light 
yield.  In  the  third  round  of  prototypes,  the  light  yield  was  restored  and  slightly  improved  by 
introducing  a  crystal  coating  layer  with  refractive  index  near  the  geometric  mean  between  that  of  the 
crystal  and  that  of  the  fiber  outer  cladding;  this  thin  layer  then  serves  as  an  antireflection  coating. 

3)  Fiber  selection  and  preparation: 

We  experimented  with  both  round  and  square  cross-section  fibers,  with  different  fluor 
formulations  and  concentrations,  and  with  diameters  ranging  from  0.5mm  to  1.0mm.  Round 
fibers  gave  nearly  a  factor  of  two  better  light  yield  relative  to  square  fibers,  which  is  due  to  much 
smaller  attenuation  losses  along  the  fibers  with  round  cross-sections.  This  arises  because  the 
square  fibers  do  not  have  a  precisely  square  cross-section,  leading  to  optical  ray  paths  which 
encounter  the  core/cladding  interface  at  more  than  the  critical  glancing  incidence  angle.  Similarly, 
1mm  fibers  have  significantly  higher  light  yield  than  0.5mm  fibers,  in  this  case  due  to  self- 
absorption  to  the  limited  Stokes  shift  between  absorption  and  emission,  combined  with  the  higher 
fluor  concentrations  required  for  efficienct  fluorescent  in  0.5mm  fibers.  We  have  used  multi-clad 
fibers  and  succeeded  in  obtaining  50%  light  yield  improvements  by  optically  coupling  aluminized 
mylar  mirrors  to  the  ends  of  polished  fiber  ribbons.  This  is  used  to  obtain  two  active  edges  and  an 
active  comer  opposite  the  readout  ends  of  the  fibers. 

41  Multilayer  stack  assembly  and  fixing: 

Multilayer  crystal/fiber  stacks  were  built  up  from  acrylic  “frames”  3.5mm  thick,  with 
cutaway  regions  on  one  side  to  guide  a  1mm  thickness  fiber  ribbon,  and  a  square  opening  in  the 
center  to  position  the  1 12mm  x  1 12mm  coated  crystal.  Fiber  ribbons  are  first  glued  to  frames,  then 
a  stack  is  built  up  by  placing  coated  crystals  onto  fibers  within  each  frame,  with  the  axes  of  the 
fiber  ribbons  alternating  between  two  orientations  as  the  stack  was  built  up.  We  devoted 
considerable  effort  to  developing  procedures  for  casting  the  modules  with  a  potting  epoxy,  so  as  to 
fill  all  interstices  without  generating  air  gaps  or  bubbles;  we  ultimately  succeeded  in  casting  nearly 
defect-free  7-crystal-deep  stacks.  Wavelength-shifting  fibers  were  also  coupled  to  the  perimeter  of 
each  crystal  in  a  stack,  with  these  fibers  read  out  separately  to  determine  the  depth-of-interaction 
within  a  stack  (these  fibers  detected  the  light  which  was  trapped  in  a  crystal  layer  and  guided  to  the 
perimeter  through  total  internal  reflection).  We  obtained  a  light  yield  of  about  10  photoelectrons 
for  0.5mm  diameter  perimeter  fibers,  which  is  more  than  sufficient  to  determine  depth-of- 
interaction  to  within  the  crystal  thickness  of  2.5mm.  The  top  of  each  module  was  epoxied  to  a 
2.5cm-thick  “light  mixer”  coupled  to  four  2”  diameter  photomultipliers  in  an  Anger  array.  The 
bottom  of  each  module  was  optically  coupled  to  reflecting  material  after  assembly,  as  were  the 
edges  of  the  light  mixer  and  those  regions  of  the  light  mixer  not  covered  by  Anger  PMTs.  Fiber 
ends  which  were  not  read  out  were  polished  and  optically  coupled  to  aluminized  mylar  mirrors. 
The  resulting  modules  have  two  “active  edges”  and  an  “active  corner”,  reflecting  the  desire  of  our 
clinical  colleagues  for  a  device  which  could  be  used  to  image  near  the  chest  wall  for  breast  imaging 
and  well  up  into  the  armpit  for  axillary  imaging. 


5)  Fiber  ribbon  routing,  multiplexing,  and  readout  coupling: 

Our  initial  tests  used  multiplexed  readout  of  two  ribbons  per  multianode  photomultiplier, 
where  a  rough  coordinate  from  the  Anger  PMT  array  was  used  to  demultiplex  between  the  two 
ribbons  on  each  fiber  readout  PMT  (i.e.  “lateral”  multiplexing  within  a  single-crystal-layer  device). 
There  is  sufficient  photocathode  area  on  our  multianode  photomultipliers  to  permit  both  4-fold 
lateral  multiplexing  and  4-fold  multiplexing  in  depth  when  using  1mm  diameter  fibers,  but  such 
heavy  multiplexing  will  not  in  future  be  needed  because  of  the  new  multianode  photomultiplier 
options  discussed  in  the  next  section.  Given  the  significant  price  of  each  multianode 
photomultiplier  it  is  important  to  keep  their  numbers  down,  but  heavy  multiplexing  decreases 
system  rate  capability  by  introducing  pile-up  and  demultiplexing  errors  at  high  event  rates.  The 
higher  rate  capability  with  new  multianode  photomultipliers  is  desirable,  but  will  require  an 
increase  in  data  acquisition  electronics  channel  count.  In  all  cases,  fibers  were  optically  coupled  to 
multianode  PMTs  by  inserting  fiber  ribbons  into  custom-built  vises,  polishing  the  fiber  ends,  then 
optically  coupling  to  the  PMTs  by  attaching  the  vises  to  custom  PMT  mounting  mechanics. 

6)  Multianode  photomultipliers  and  associated  electronics: 

We  used  Hamamatsu  R5900-L16  multianode  photomultipliers  to  read  out  the  fibers,  where 
each  multianode  PMT  contained  16  separate  photocathode  regions  each  16mm  x  1mm.  The  16 
anodes  were  resistively  coupled  in  a  chain  and  read  out  through  charge  division.  We  were  able  to 
obtain  spatial  resolution  of  less  than  1mm  FWHM  for  8  photoelectrons  incident  through  a  single 
1mm  fiber;  our  spatial  resolution  is  therefore  not  limited  by  our  PMT  readout,  but  rather  by 
upstream  optics.  The  two  charge-division  outputs  (whose  ratio  give  the  transverse  coordinate  for 
light  incident  through  a  fiber  ribbon,  and  whose  sum  gives  the  total  fiber  light  yield)  were  read  out 
through  custom-built  4-PMT  bases  and  coupled  to  xlO  preamplifiers  which  we  also  constructed. 
The  single  photoelectron  peak  from  the  multianode  photomultipliers  was  very  pronounced, 
simplifying  detector  gain  calibration  considerably. 

7)  Coincidence  trigger: 

We  used  commercial  NIM  electronics  to  construct  a  simple  but  adequate  trigger  for  our 
system.  It  begins  with  an  analog  fan-in/fan-out  module  which  combines  the  four  Anger  PMT 
signals  (which  are  split  in  two  by  50-ohm  splitters,  with  the  second  half  going  to  ADCs)  to  make  a 
single  signal  from  each  module;  these  signals  are  then  discriminated,  and  a  trigger  resulting  from  a 
coincidence  between  the  two  discriminated  signals  generates  a  gate  which  initiates  ADC  conversion 
and  later  event  readout.  We  have  assembled  a  second  commercial  NIM  trigger  subsystem  as  part 
of  the  portable  system  we  are  building. 

8)  Data  acquisition  electronics  and  computer  interface: 

We  used  commercial  LeCroy  2249W  ADCs  to  digitize  the  outputs  of  both  4  Anger  PMTs 
per  module  and  16  multianode  PMT  charge-division  outputs  per  module.  A  trigger  was  generated 
from  the  coincidence  of  discriminated  output  of  analog-summed  Anger  PMT  signals,  and  each  of 
the  40  digitized  signals  was  delayed  by  200ns  to  land  within  a  lOOOns-width  gate.  A  second, 
lower  threshold  was  also  set  on  each  Anger  PMT  sum,  and  the  outputs  of  these  discriminators  sent 
to  a  LeCroy  2228  TDC  to  determine  the  coincidence  time  resolution  as  a  function  of  threshold. 
CAMAC  readout  of  the  42  detector  channels  was  performed  with  a  crate  controller  designed  at  the 
Budker  Institute  for  Nuclear  Physics  in  Siberia,  which  is  capable  of  nearly  lOKHz  readout  of  this 
many  channels  and  is  coupled  to  a  PC.  In  addition  to  the  2249W  ADC  system  we  have  tested  a  few 
modules  of  FERA-Fast  Encoding  and  Readout  ADCs.  These  zero-suppressing  and  fast 
conversion-time  modules  should  extend  our  system  rate  capability  to  50-100  Khz,  which  is  a 
higher  rate  than  we  expect  at  reasonable  patient  doses. 


Figure  13:  Drawing  #10  —  Conceptual  drawing  of  two  proposed  detector  modules  mounted  on 

mammograph  ~  note  that  actual  mammograph  (unlike  unit  shown) 
has  a  detachable  bucky/cartridge  support  structure. 
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Abstract 


We  have  characterized  the  new  Hamamatsu  R5900-L16 
multianode  photomultiplier  with  regard  to  its  suitability  for 
position-sensitive  readout  with  just  two  analog  channels  for 
each  16-anode  PMT.  Because  of  the  low  capacitance  of  the 
anode  structure  of  the  PMT,  we  obtain  an  output  pulse  width 
of  less  than  20  ns  for  narrow  input  pulses.  We  were  able  to 
obtain  position-sensitive  readout  of  single  photoelectron 
pulses  with  a  spatial  width  of  less  than  1  mm  FWHM  for  a 
collimated  light  source  and  about  2  mm  FWHM  for  single 
fiber.  Combined  with  the  good  single  photoelectron  peak  for 
R5900-L16,  this  provides  excellent  performance  at  very  low 
light  levels.  We  discuss  application  of  this  device  to  the 
precise  determination  of  gamma-ray  interaction  position 
within  thin  crystal  scintillators. 

I.  Introduction 

Position  sensitive  photomultipliers  have  been  used  by  a 
number  of  groups  for  the  readout  of  scintillating  or 
wavelength-shifting  fibers.  Crossed-wire  and  multianode 
PMTs  have  also  been  used  for  the  position-sensitive  readout  of 
extended  light  sources  [2-5].  At  low  light  levels,  the  cross¬ 
talks  within  the  PMT  can  result  in  significant  degradation  in 
spatial  resolution.  There  are  several  sources  of  such  cross¬ 
talks:  light  spreading  through  the  entry  window  glass  on  the 
way  to  the  photocathode,  imperfect  electron  optics  between 
photocathode  and  dynode  structure,  secondary  electron  shower 
spreading  within  the  multiplier  assembly,  or  just  electrical 
cross-talks  between  anodes  (induced  signals).  A  significant 
difference  between  position-sensitive  cross-wired  PMTs  like 
the  Hamamatsu  R2486,  versus  multianode  PMTs  like  the 
Hamamatsu  R4760,  is  the  considerably  higher  cross-talks  and 
therefore  decreased  spatial  resolution  at  low  light  level  of  the 
former.  The  Hamamatsu  R2486  for  instance  requires  10000 
photons  at  490  nm  (about  1000  photoelectrons)  to  obtain  a 
spatial  resolution  of  0.2  mm  FWHM  [2];  the  resolution  is 
proportional  to  the  square  root  of  the  number  of  photoelectron 
per  event.  Another  difficulty  with  cross-wired  PMTs  for  some 
applications  is  its  limited  rate  capability  due  to  significant 
pulse  width  with  typical  resistor  chain  readout. 
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Nonetheless,  the  reduced  electronics  channel  count  and 
lower  device  cost  (particularly  per  unit  of  useful  photocathode 
area)  associated  with  crossed-wire  PMTs  have  made  them 
attractive  for  several  applications,  including  for  the  readout  of 
sets  of  scintillating  fibers. 

There  is  some  success  in  attempts  to  use  charge  division 
method  for  two-dimensional  readout  of  an  array  of  optical 
fibers  by  multianode  PMT  having  just  four  analog  channels 
[5].  But  again,  one  needs  hundreds  of  photons  to  obtain 
submillimeter  resolution. 

We  investigated  the  properties  of  the  newly  developed 
PMT  -  the  Hamamatsu  R5900-L16  -  because  of  its  very 
promising  specifications  caused  mainly  by  recently  developed 
microchannel  dynode  structure. 

II.  Experimental  Setup 

A.  R5900-L16  -  multianode  PMT 

The  photomultiplier  described  here  has  external 
dimensions  of  merely  a  cube  25*25*25  mm.  It  has  flat 
bialkali  photocathode,  10-stage  dynode  structure  and  16  anode 
stripes  0.8  *  16  mm  regularly  spaced  at  1  mm  pitch.  An 
important  feature  of  the  PMT  distinguishing  it  from  previous 
ones  is  that  the  microchannel  dynodes  work  effectively  like 
separate  PMTs  in  terms  of  electron  optics.  This  minimizes 
electron  shower  smearing  and  allows  one  to  obtain  very  fine 
spatial  resolution  as  will  be  shown  below.  Another  useful 
feature  of  the  microchannel  structure  is  a  relatively  high  gain 
at  low  operating  voltage  -  about  3..4*106  at  800V  .  The  first 
prototype  of  R5900-L16  sent  to  us  by  Hamamatsu  had  a  size 
of  photocathode  of  15.5  *  15.5  mm  which  is  a  0.5mm 
narrower  than  dynode  structure.  Most  of  the  measurements 
were  done  on  this  particular  PMT.  The  photocathode  size  of 
the  last  version  of  R5900-L16  is  significantly  larger  - 
approximately  21  *  22  mm.  We  have  repeated  some  of  the 
measurements  to  check  the  impact  of  this  design  change  on  the 
performance. 

B.  Charge  Division  Readout 

We  have  connected  16  anodes  of  the  PMT  by  a  linear  array 
of  15  100-ohm  resistors  to  couple  16  individual  outputs  into  2 
charge  division  outputs.  Each  of  the  two  signals  was  directly 
coupled  through  50-ohm  cables  to  a  LeCroy  2249W  operating 
with  a  150  ns  gate  width.  To  minimize  artifacts  caused  by 
discrete  nature  of  ADC  data  we  have  inserted  10-fold 
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amplification  of  the  charge  division  signals  giving  us  a  single 
photoelectron  peak  at  about  35  channels  of  the  ADC  over  the 
pedestal. 

B.  Data  Acquisition  System 

For  measurements  described  here  we  used  rather  simple 
data  acquisition  system  consisting  of  a  CAMAC  ADC 
mentioned  above  and  micro-VAX  computer  under  Open  VMS 
connected  to  the  CAMAC  crate  and  controller  via  SCSI-bus. 
Data  collected  on  a  disk  of  the  computer  then  were  analyzed  in 
a  PAW  (CERN  written  data  analysis  software)  session. 

III.  Tests  and  Measurements 

A.  Timing  and  Single  Photoelectron 
Characteristics 

The  PMT  itself  (when  its  anodes  are  read  out  separately) 
has  fast  rise  time  and  narrow  pulse  width  -  0.6  ns  rise  time.  1 .7 
ns  fall  time  when  loaded  to  50  Ohm.  After  connecting  all 
anodes  by  a  linear  array  of  100-Ohm  resistors,  we  obtain  a 
pulse  width  of  less  than  20  ns.  To  check  single  photoelectron 
amplitude  peak  we  connected  together  two  ends  of  the  charge 
division  chain.  The  PMT  was  illuminated  by  a  green  LED 
[  giving  0. 1  photoelectron  in  average.  The  distribution  is  shown 
*  on  Fig.  1.  The  peak-to- valley  ratio  of  almost  3:1  is  a  very 
-  desirable  characteristic  for  low  light  level  applications. 
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Figure  1:  Single  phoioelectron  peak  of  the  R5900-L16.  All  anodes 
are  read  out  together.  HV=800V. 

B.  Charge  Division  Calibration 

Ideally,  if  we  know  parameters  of  the  chain  -  we  can 
calculate  the  calibration  coefficient  between  charge  division 
ratio  R=(A2-A1)/(A1+A1)(A1  and  A2  -  amplitudes  measured 
from  two  ends  of  the  chain)  and  position  in  the  chain, 
considering  each  anode  to  be  a  current  source.  The 
correspondence  between  charge  division  node  number  and 
measured  ratio  R  should  be  linear.  Therefore  even  if  several 
photoelectrons  are  distributed  over  several  PMT  anodes  the 
measured  ratio  R  would  correspond  to  the  center-of-gravity  of 


received  photoelectrons.  In  our  particular  case  we  have 
Rin=50  Ohm  and  15  resistor  of  the  chain  100  Ohm  each.  So 
ideally  for  point  source  from  the  first  or  the  last  anode  we 
should  measure  the  ratio 

R=+-(50+ 15*1 00-50)1(50+ 15*1 00+50)=+-0.9375 . 

In  real  situation  one  has  to  correct  this  ratio  for  different 
gains  and  input  impedances  of  the  two  amplifiers.  Practically 
we  found  that  the  ratio  R  ranges  from  -1.05  to  1.05.  which 
means  that  anodes  others  than  illuminated  by  light  source 
produce  induced  signals  of  inverse  polarity.  Indeed  we  have 
seen  and  measured  such  signals  on  the  scope  -  the  sum  of  all 
anodes  but  illuminated  one  gives  approximately  15%  of  the 
main  signal.  Neglecting  the  fact  of  induced  signals  we  could 
calculate  the  linear  coefficient  for  conversion  of  the  measured 
ratio  R  to  the  signal  center-of-gravity  position: 

K=  1 5/(  1 .05-(- 1 .05))=7 . 14, 

where  15mm  -  separation  between  the  first  and  the  last 
anodes.  Using  this  coefficient  we  have  plotted  then  center-of- 
gravity  coordinate  versus  actual  position  of  the  light  source. 


actual  position,  mm 


Figure  2:  Reconstructed  coordinate  lineanty.  A  0.5  mm  fiber 
coupled  to  a  green  LED  was  positioned  by  a  calibrated  translauon 
stage  with  an  increment  of  0.5  mm.  At  each  position  corresponding 
charge  division  coordinate  was  calculated  and  plotted.  Light  level  is 
approximately  100  photoelectrons. 


C.  Intrinsic  Spatial  Resolution  of  the  PMT 

In  order  to  measure  intrinsic  resolution  of  the  PMT  we 
have  used  collimated  parallel  light  beam  produced  by  a  green 
LED  and  a  black  shield  with  two  holes  of  0.3mm  diameter. 
The  shield  was  mounted  directly  on  the  PMT  window  so  we 
were  able  to  avoid  light  spreading  in  the  window  glass  which 
thickness  is  about  1 .5  mm.  The  LED  intensity  was  adjusted  to 
give  about  0.1  photoelectrons  per  pulse.  Only  events 
corresponding  to  single  photoelectron  peak  were  recorded 
therefore  converted  photon  should  have  coordinate  of  the  first 
or  the  second  hole  of  the  collimator.  The  distribution  and  two- 
gaussian  fit  are  shown  on  Fig.3.  The  reconstructed  distance 
between  two  peaks  well  matches  the  hole  separation  of  1mm. 
The  width  of  each  of  the  two  peaks  is  about  0.25mm  RMS  (or 
0.6mm  FWHM). 
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Figure  3:  Reconstructed  coordinate  distribution  for  events  selected 
around  single-photoelectron  peak  +/-  1  RMS.  Light  source  is  a 
double  beam  produced  by  green  LED  through  a  collimator  -  black 
film  with  two  holes  of  0.3  mm  diameter  separated  by  1  mm.  Two 
peaks  clearly  separated  at  a  distance  of  about  1mm. 

D.  Single  Fiber  Response 

There  are  two  difference  between  collimated  beam  and 
light  coming  from  a  fiber.  First,  fiber  photons  leave  the 
surface  of  the  fiber  uniformly  over  the  its  cross-section. 
Second,  these  photons  have  an  angular  spreading  of  about  20- 
30  degrees,  depending  on  type  of  the  core  and  the  cladding 
materials.  These  two  factors  together  broaden  the  distribution 
of  the  photons  converted  on  the  photocathode.  In  order  to 
evaluate  this  effect  we  have  repeated  previous  measurement 
using  instead  of  the  collimated  light  beam  -  0.5mm  clear  fiber 
coupled  to  a  green  LED  and  to  the  PMT  by  optical  grease. 
The  distribution  is  shown  on  Fig.4.  It  doesn't  much  differ  from 
the  distribution  on  Fig.3,  which  means  that  photoelctrons  in 
this  particular  test  are  mostly  distributed  over  two  anodes  or 
within  2  mm  in  physical  space.  Another  useful  test  is  how  the 
resolution  depends  on  number  of  registered  photoelectrons. 
We  found  that  the  spatial  resolution  scales  with  the  inverse  of 
the  square  root  of  the  number  of  the  photoelectrons  to  well  in 
excess  of  20  photoelectrons.  The  results  for  different  light 
levels  are  shown  on  Fig.5. 
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Figure  4:  Reconstructed  coordinate  distribution  for  all  events  with 
amplitude  greater  than  0.5  photoelectron..  Light  source  is  a  green 
LED  coupled  to  the  PMT  through  a  0.5mm  clear  fiber  giving  an 
average  of  0.1  photoelectron. 


Figure  6:  Visible  light  profile  across  the  photocathode  of  the  PMT. 
One  end  of  a  0.5  mm  fiber  is  coupled  to  a  green  LED,  another  end 
was  moving  in  a  proximity  of  the  PMT  face  with  a  calibrated 
translation  stage. 

E.  Field  of  view  test . 

All  previous  results  were  obtained  in  just  few  positions  on 
the  PMT  window.  Important  questions  also  are  how  uniform 
and  how  big  is  field  of  view  (photocathode  area).  Some  edge 
effects  are  visible  on  Fig.2  so  useful  area  for  the  first  PMT 
tested  was  about  14*  14mm.  The  result  of  a  uniformity  test  is 
shown  on  Fig.6.  The  width  of  "uniform"  area  is  about  15mm 
if  one  would  agree  to  lose  50%  of  the  signal.  This  lost  would 
mean  less  number  of  photoelectrons  registered  and 
consequently  some  degradation  in  resolution  on  the  edge. 

As  it  was  mentioned  before,  the  last  version  of  the  PMT 
has  wider  active  photocathode  area.  Therefore  we  have 
repeated  the  measurements  reflected  on  Figs.2  and  6  and 
found  that  "useful"  area  of  the  last  version  of  the  PMT  (based 
on  the  same  criteria)  is  about  1mm  wider  in  both  directions 
(i.e.  16*  16mm). 


IV.  Discussion 

We  are  particularly  interested  in  application  of  this  PMT  to 
the  readout  of  wavelength-shifting  fibers  coupled  to  thin 
layers  of  scintillator  crystals,  for  use  in  a  PET  detector  were 
designing  [6].  In  this  application,  the  spread  of  the  light  across 
a  ribbon  of  photomultiplier  has  a  FWHM  about  equal  to  the 
thickness  of  the  crystal  -  2  mm  for  our  device.  Achievable 
spatial  resolution  limited  by  this  spread  could  be  in  principle  a 
FWHM  of  2  mm  divided  by  the  square  root  of  the  number  of 
photoelectrons  measured.  This  resolution  is  then  further 
degraded  by  the  spatial  resolution  of  the  fiber  readout  PMT, 
which  from  the  above  measurements  is  about  twice  better. 
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Figure  5:  Reconstructed  coordinate  distribution  for  different  light 
levels  (upper  row  from  left  to  the  right:  I,  2  and  3  photoelectrons; 
bottom  row  from  left  to  the  right:  5,  10  and  20  photoelectrons).  Light 
is  coming  from  a  green  LED  coupled  to  the  PMT  through  a  0.5mm 
clear  fiber. 

Since  the  two  effects  are  incoherent,  we  could  preserve  an 
ntnnsic  resolution  of  our  crystal-fiber  readout  design 
mdisturbed  by  the  PMT  itself.  At  a  current  cost  of  less  than 
$1500  in  single-unit  quantities,  this  PMT  is  cost-effective  in 
comparison  with  competing  devices,  particularly  in 
multiplexed  geometries  such  as  one  obtains  when  reading  out 
me  end  each  of  many  fibers  at  each  PMT  anode  input  (e.g.  14 
ibbons  of  lmm-diameter  fibers  each  14mm  wide,  for  total  of 
196  fibers,  may  be  read  out  with  just  two  PMTs  by  differing 
:he  orientation  of  the  readout  PMTs  by  90  degrees  about  the 
fiber  axis  at  the  two  fiber  ends  -  presuming  that  just  1  fiber  of 
the  196  fibers  on  a  given  event).  The  level  of  multiplexing 
possible  at  low  light  levels  is  critically  dependent  on  the 
spatial  resolution  achievable. 


v.  Conclusion 

We  have  studied  the  properties  of  a  photomultiplier  with 
charge  division  readout  suitable  for  scintillating  or 
wavelength-shifting  fiber  applications  with  submillimeter 
resoluuon  at  low  light  levels  of  just  few  photoelectrons.  It  was 
shown  that  intrinsic  resolution  of  the  PMT  is  better  than  1mm 
FWHM  even  for  single  photoelectron  signals.  The  spreading 
of  the  light  coming  from  the  fibers  to  the  PMT  windows 
deteriorates  resolution  greater  than  electron  optics  does.  From 
the  other  hand  we  have  found  that  electrical  cross-talks 
between  anodes  don’t  disturb  achievable  resolution  and 
linearity  of  the  charge  division  method. 
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Abstract 

A  low-cost  replacement  for  popular  LeCroy  2249  series 
charge-integrating  ADCs  with  more  than  100  times  faster  data 
throughput  is  designed  and  tested.  The  A/D  converter  chips 
used  in  the  design  provide  700  nsec  conversion  time  and  25 
MHz/32-bit  parallel  bus  made  of  a  twisted  pair  ribbon  cable 
connected  to  a  high-speed  buffered  PCI  interface  is  capable  of 
100  MByte/sec  data  transfer  to  an  ordinary  IBM  PC.  The 
front-end  of  each  ADC  channel  is  equipped  with  an  unipolar 
Baseline  Restoration  Circuit  allowing  use  of  AC  coupling  at 
very  high  signal  rates  (RF  preamplifiers,  HV  decoupled  signal 
sources  like  PMTs,  wire  chambers,  etc.). 

Being  developed  for  Positron  Emission  Tomography,  the 
system  could  also  be  used  for  a  range  of  High  Energy  Physics 
applications.  Spectroscopy,  Medical  Imaging  and  others  where 
high  event  rate  occurs. 

I.  INTRODUCTION 

This  paper  presents  the  design  and  some  performance 
results  for  the  newly  developed  ADC  sytem  for  use  in  our 
Positron  Emission  Tomography  Project.  The  need  for 
developing  such  a  system  arises  from  various  limitations  of 
commercially  available  A/D  converters.  For  this  particular 
project  we  wanted  a  portable  system  of  about  a  hundred  ADC 
channels  with  a  moderate  10-bit  resolution  but  capable  of 
digitizing  and  transferring  the  data  at  up  to  50k  events  per 
second  rates.  Systems  satisfying  our  requirements  exist  but 
they  employ  very  expensive  and  bulky  standards  like 
FASTBUS.  Other  less  expensive  solutions  including  LeCroy 
FERA[1]  system  are  acceptable  in  terms  of  conversion  rate 
and  resolution  but  the  speed  of  data  transfer  is  too  low  for  our 
application.  Handling  an  event  rate  much  higher  than  few  kHz 
represents  a  serious  technical  challenge.  For  instance,  in  order 
to  transfer  200  byte  events  at  50  kHz  rate  one  needs  a 
computer  interface  with  at  least  10  Mbytes/sec  data 
throughput.  Fortunately,  recent  progress  in  CMOS  technology 
and  Programmable  Logic  Devices  (PLDs)  allows  one  to 
design  hardware  meeting  the  above  requirement  at  a 
reasonable  cost.  Also  a  range  of  A/D  converter  chips  with  10- 
12-bit  resolution  and  conversion  time  shorter  than  1  usee  is 
now  available  from  different  vendors  for  $8-15.  In  addition  to 
this,  In-System  Reprogrammability  of  last-generation  PLDs 
makes  the  debugging  stage  of  the  design  cycle  much  shorter  in 
time  and  less  labor-intensive.  All  these  reasons  led  us  to 
developing  our  own  ADCs  with  a  PC-compatible  interface. 


H.  System  Architecture 

A.  General  Description 

We  have  chosen  to  implement  the  first  version  of  the 
design  in  the  CAMAC  standard  although  we  don't  use  the 
CAMAC  backplane  other  than  for  the  power  supply.  The  idea 
was  to  be  able  to  switch  to  a  CAMAC-like  crate  containing 
just  a  moderate-wattage  power  supply  in  it,  and  get  rid  of  the 
bulky,  heavy  and  expensive  full-featured  CAMAC  crate  when 
the  system  has  to  be  transported  between  hospitals  for  clinical 
trials.  For  the  same  reason  it  was  convenient  to  have  a  trigger 
module  residing  in  the  crate  right  next  to  the  ADC  cards  and 
have  no  extra  NIM  crate/modules. 

In  order  to  get  the  ADC  data  to  the  computer,  we  used  a 
buffered  PCI  interface  compatible  with  an  ordinary  PC,  which 
we  had  previously  developed  for  another  project.  More 
detailed  descriptions  of  the  trigger  unit  and  the  interface  will 
follow  in  two  separate  paragraphs  below. 

B.  ADC  card 

We  decided  to  lay  out  10  ADC  channels  with  data 
interface  per  single  CAMAC  card  although  we  had  probably 
enough  space  to  accommodate  16  channels.  This  choice  was 
dictated  by  the  number  of  I/O  pins  of  the  chosen  PLDs.  For 
our  implementaton  it's  convenient  to  have  delay  lines 
embedded  in  the  design,  so  the  user  doesn’t  have  to  externally 
delay  analog  signals.  We  used  200  nsec  solid  state  delay  lines 
from  JBM[2]  to  equip  each  ADC  channel. 

The  12-bit  data  from  all  10  channels  are  latched  in  parallel 
in  two  Altera  [3]  EPM7128S  PLDs  to  be  then  transferred  to 
the  computer  host.  In  order  to  achieve  maximum  possible  data 
transfer  rate  to  a  computer  and  yet  to  keep  the  system 
reconfigurable  we  have  chosen  a  unidirectional  daisy-chain 
priority  arbitration  between  ADC  cards.  It  works  as  following: 
the  first  card  has  the  highest  priority,  and  as  soon  as  it  finishes 
its  data  transfer  through  the  parallel  32-bit  bus,  the  second 
card  gains  access  to  the  bus  and  sends  its  data.  Then  the  third 
card  starts,  and  so  on.  The  same  PLDs  mentioned  above 
provide  all  control  signals,  daisy-chain  arbitration  and  self¬ 
locking  of  the  system  (to  prevent  signal  jamming  during  A/D 
conversion). 

C.  More  Detailed  Description  of  one  ADC  channel 

One  serious  problem  with  a  high  signal  rate  is  the  so- 
called  baseline  shift,  or  floating  of  DC  offsets  when  the  signal 
rate  fluctuates.  This  has  to  be  eliminated  before  the  signal 
integration  so  that  the  integral  only  represents  pulse  area.  In 
order  to  design  an  integrator  suitable  for  accurately  integrating 


small,  fast  PMT  pulses  which  are  superimposed  on  a  DC 
offset  voltage  from  decoupling  capacitors,  several  issues  are  to 
be  addressed.  The  first,  is  to  have  sufficient  signal 
preamplification  so  that  the  integrating  circuit  can  develop  a 
reasonable  output  voltage  for  A/D  conversion.  Second,  we 
have  to  compensate  for  the  effect  of  changing  DC  offset  on 
the  integrated  output,  or,  in  other  words,  provide  baseline 
restoration  (BLR).  And  third,  we  have  to  use  a  sufficiently  fast 
op-amp  in  the  integrator  that  non-linear  effects  are  not 
significant.  Also  the  op-amp  has  to  have  very  low  input 
current  in  order  to  minimize  drooping  of  the  integrated  signal. 

Signal  preamplification  requires  the  use  of  very  wideband 
amplifiers  if  signal  fidelity  is  to  be  preserved.  Although  such 
amplifiers  exist,  they  tend  to  exacerbate  the  DC  offset  problem 
because  they  usually  have  larger  voltage  offsets  and  input  bias 
currents.  A  better  approach  is  to  use  shaping  amplifiers  rather 
than  linear  ones,  thereby  allowing  the  use  of  lower  bandwidth 
(and  lower  cost)  amplifiers.  If  the  shaping  uses  a  single 
identical  lag  for  each  amplifier,  then  the  input  pulses  will  be 
stretched  without  adding  overshoot.  Lengthening  the  fast  input 
pulses  also  allows  us  to  use  a  slower  op-amp  for  the 
integrator. 

Baseline  restoraton  before  the  integrator  input  is  relatively 
simple  if,  as  in  our  case,  the  pulses  to  be  integrated  are 
unipolar.  The  circuit  that  removes  the  offset  is  usually  called  a 
baseline  restorer  and  can  be  either  passive  or  active.  For  the 
lowest  offset,  an  active  circuit  is  preferred.  In  our 
implementation,  an  OPA660  (Burr-Brown)  transconductance 
amplifier  senses  the  offset  voltage  to  develop  a  feedback 
current  that  biases  the  preamplifier  input  in  a  way  that  the 
output  offset  is  reduced  to  nearly  zero.  A  diode  is  used  in 
series  with  the  output  of  the  transconductance  amplifier  so  that 
only  one  polarity  of  offset  is  corrected.  The  very  high 
bandwidth  of  OPA660  (850  MHz)  is  required  so  that  the 
baseline  restoration  acts  only  on  the  instantaneous  offset  rather 
than  on  the  average.  The  final  version  of  our  BLR  circuit 
shows  a  very  small  offset  (less  than  1%  of  the  input 
amplitude)  when  an  AC  coupled  pulser  feeds  the  ADC  input 
with  a  rectangular  signals  of  50  nsec  width  at  a  rate  of  5  MHz 
(i.e.  signal  is  present  25%  of  the  time). 

The  op-amp  used  as  the  feedback  integrator  (we  used 
LT1361  with  JFET  inputs  from  Linear  Technology[3])  must 
have  wide  bandwidth  so  that  its  input  summing  junction 
remains  at  zero  volts  at  all  times  as  was  mentioned  above.  If 
the  bandwidth  is  insufficient,  then  the  summing  junction  will 
move  with  the  input  signal  and  the  integral  will  no  longer  be  a 
linear  function  of  the  input  pulse  area.  By  using  pulse  shaping 
in  the  input  amplifier,  the  required  bandwidth  of  the  integrator 
is  reduced,  allowing  a  wider  choice  in  op-amps.  Except  during 
the  integration  interval,  typically  of  order  of  hundreds 
nanoseconds,  the  integrator  output  is  held  at  zero  volts  by  two 
DMOS  switches  (SST215),  one  shorting  the  integrating 
capacitor  and  the  other  shorting  the  input  to  the  integrator. 
These  are  extremely  fast  switches  (one  nanosecond  on/off 
transition)  which  have  to  be  are  turned  off  to  allow 
integration. 


For  the  A/D  converter  itself  we  decided  to  use  LTC1410 
(also  from  Linear  Technology,  see  [5])  -  capable  of  1.25M 
conversions  per  second  with  12  bit  resolution.  We  don't 
however  use  the  power  saving  and  microprocessor-type  data 
output  control  features  offered  by  the  manufacturer. 

D.  Computer  Interface 

Technical  details  of  the  interface  design  lie  definitely 
beyond  the  scope  of  this  paper  but  a  general  description  will 
be  useful.  The  interface  is  a  single  PCI  board  residing  in  an 
IBM  PC  system  block  (currently  the  computer  runs  the  Linux 
operating  system  but  we  are  considering  to  write  a  driver  for 
the  interface  under  Windows  95/98  OS).  It  accepts  32  bits  of 
data  (TTL  levels)  in  parallel,  accompanied  with  a  strobe 
signal,  all  of  the  above  transmitted  through  a  twisted  pair 
ribbon  cable  of  up  to  10  m  in  length.  Sequential  data  words 
have  to  have  a  minimum  time  separation  of  20  nsec  between 
them,  providing  a  peak  throughput  of  the  interface  itself  of 
200  Mbytes/sec.  In  our  particular  implementation  we’ve 
chosen  40  nsec  separation  of  sequential  words  within  one 
event  which  saturates  the  data  bus  at  100  Mbytes/sec. 
Certainly  such  an  enormous  data  flow  can  neither  be  recorded 
to  a  disk,  nor  be  processed  on  the  fly  (at  least  with  a  single 
processor).  Our  target  data  flow  was  rather  several  Mbytes/sec 
and  can  be  handled  by  a  moderate  IBM  PC.  In  order  to 
eliminate  dead  time  when  the  CPU  reads  the  data  out,  the  PCI 
interface  has  two  memory  banks:  when  the  first  is  full  and  is 
being  read  by  the  computer,  the  second  bank  is  available  for 
writing  into  it  from  the  ADCs. 

In  the  case  of  our  initial  implementation  of  5  10-channel 
cards,  the  whole  data  transfer  lasts: 

5  x  10  x  Vi  x  40  nsec  =  1  microsec, 

(  Vi  is  here  because  data  words  get  transfered  by  pairs). 
Combined  with  the  conversion  time  of  700  nsec  and  small 
daisy-chain  overhead  it  yields  a  maximum  throughput  of 
about  500k  events  per  second  if  the  computer  is  powerful 
enough.  For  a  greater  number  of  ADC  cards  the  maximum 
throughput  will  be  decreased  linearly. 

E.  Trigger  Board/Gate  Generator 

Strictly  speaking,  the  above  ADC  cards  require  only  one 
external  signal  -  TTL-leveled  “Gate”  determining  how  long 
input  analog  signals  are  to  be  integrated.  The  end  of  “Gate” 
signal  starts  A/D  conversion  in  all  ADC  cards  simultaneously. 

We  found  it  convenient  to  combine  discriminators, 
coincidence  circuit  and  gate  generator  (with  apropriate  gate 
delay  and  width  adjustments)  in  one  CAMAC  card,  together 
with  TTL  pulser  for  calibration  needs.  All  logic  of  the  above 
functions  fits  within  just  one  inexpensive  Altera  PLD  chip  - 
EPM7064SQC44.  Currently  we  have  four  independent 
discriminator  channels  (based  on  LeCroy  discriminator  chips  - 
MVL107S  [1],  with  ouputs  converted  from  ECL  to  TTL). 
They  work  reliably  with  a  minimum  threshold  of  about 
20  mV.  The  width  and  the  delay  of  the  “Gate”  signal  are 
easily  programmable  in  order  to  satisfy  any  particular 
requirements. 


Figure  I .  Example  of  amplitude  distributions  for  four  randomly  chosen  channels.  The  pedestal  spike  widths  are  about  1.3  LSBs  FWHM. 
Acquired  with  R5900  FMTs  and  lOx  preamplifiers  attached  to  the  above  ADCs.  Average  signal  is  very  small  (~0.8  photoelectrons)  so  single¬ 
photoelectron  peaks  are  clearly  visible. 


Again,  because  of  the  In-System  Reprogrammability 
feature  of  these  Altera  PLDs,  it’s  easy  to  reconfigure  the  logic 
of  the  trigger,  say,  for  instance  from  total  OR  of  all  input 
signals  to  some  more  complicated  coincidence  of  their  pairs, 
etc. 

m.  Performance  Tests  and  Results 

The  very  first  check  we  have  undertaken  after  initial 
debugging  of  the  design  was  to  see  if  the  data  gets  corrupted 
during  transfer  from  ADCs  to  the  computer.  Fourtunately  we 
didn’t  find  any  damaged  datawords  among  many  millions  of 
recorded  bytes.  It  proves  that  the  daisy  chain  and  twisted  pair 
databus  are  functioning  properly. 

A  second  important  test  is  determining  the  maximum 
speed  of  the  system.  We  found  that  a  single  processor 
Pentium-II/266  MHz  PC  is  capable  of  recording  of  about  35k 
events  (100-byte  size)  per  second  to  its  hard  disk.  This  is  a 
pure  disk  limitation,  since  the  speed  of  just  reading  the  data 
into  PC  memory  exceeded  100k  events  per  second. 

We  measured  the  widths  of  pedestal  distributions  in 
presence  of  high  voltage  on  PMTs  with  lOx  preamplifiers 
connected  to  several  ADC  channels  (see  Fig.l).  The  average 
width  was  about  1.3  LSBs  FWHM. 

We  didn’t  spend  a  lot  of  time  characterizing  overall  non- 
linearities,  but  we  did  check  the  differential  non-linearity  in 
several  channels  by  acquiring  very  high  statistics  of  smooth 
amplitude  distribution.  In  good  agreement  with  the 
manufacturer’s  specification  [5],  non-smoothness  of  measured 
distributions  did  not  exceed  30%. 

Also  we  measured  long  and  short-term  drifts  of  the  ADCs 
and  found  that  the  pedestals  in  all  channels  shifted  by  1-2 
LSBs  (Least  Significant  Bits)  within  1  minute  after  turning  on 
the  power  supply.  After  this  first  minute,  the  values  in  most  of 
the  channels  stayed  stable  vs  time  and  temperature  (we  have 
repeated  the  measurements  for  several  weeks)  within  +/-1 
LSB. 

Among  other  parameters  characterizing  ADC  perfomance 
are: 

integration  gate  width  -  20  nsec  -  20  usee; 
sensitivity  -  0.2  pC/LSB; 
input  impedance  -  50  Ohm. 


Power  requirements  per  ADC  board: 

+6V  -  0.5 A,  -6V  -  0.4A,  +/-  12V  -  0.12A,  or  about  8.5W  of 
total  power  consumption.  As  one  can  see,  the  heat  dissipation 
is  not  substantial  so  it  doesn’t  need  forced  airflow  for  the 
cooling. 
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Absiraa 

Bv  usine  wavelength-shifting  libers  coupled  to  thin  plates 
of  LSO  and  to  photomultipliers,  we  have  demonstrated  spatial 
resolution  of  2mm  FWHM  lor  photocapturc  events. 
Optimized  systems  are  believed  to  be  capable  ot  better  than 
1mm  FWHM  resolution,  without  requiring  cither  scry  small 
crystals  or  large  numbers  of  photosensors.  With  this  sensor 
resolution,  annihilation  acollincaritv  dominates  system 
resolution  tor  large  PET  rings.  By  radially  stacking 
alternating  layers  of  thin  crystals  and  wavelength-shifting 
fibers,  a  direct  depth-of-interacnon  measurement  is  possible. 
Energy  measurements  and  coincidence  timing  may  he 
provided  bv  intzger  photomultipliers  in  a  dual-photodciccior 
geometry  Results  from  computer  simulations  and  component 
response  measurements  on  proot-ot*conccpt  prototypes  are 
presented. 

I.  INTRODUCTION 

We  are  developing  very  high-precision  PET  detector 
modules  with  depth-ol-intcraction  sensitmts.  using  a  new 
design  incorporating  wavelength-shifting  liber  readout  ol  thin 
LSO  crystals.  By  building  modules  from  thin  ,r\stal  layers 
stacked  radiallv.  one  can  obtain  both  depth-oi-mteraction 
information  and  line  spatial  resolution  in  the  tangential  and 
axial  directions.  Component-level  testing  ot  this  technique 
has  demonstrated  <2  mm  FWHM  spatial  resolution  lor  2mm 
thick  LSO  crystals.  Optimization  and  extrapolation  to  full 
PET  systems  may  yield  devices  with  resolution  as  small  as 
1mm  FWHM  using  1 -2mm  thick  crystal  layers:  the  resolution 
is  expected  to  vary  linearly  with  the  crystal  lavcr  thickness. 
Depth-ot-interaction  sensitivity  is  essential  lor  uniformly 
high-precision  PET  throughout  a  large  ticld-ot-view.  it  one  is 
to  preserve  the  detector  thickness  needed  tor  high  efficiency 
and  high  event  statistics.  Dcpth-ot -interaction  sensitivity  also 
minimizes  the  necessary  detector  diameter,  thcrctiv  decreasing 
both  acollincantv  errors  and  system  costs,  while  increasing 
solid  angle  coverage.  Our  modular  device  design  permits 
testing  of  a  variety  of  detection  geometries,  including 
geometries  appropriate  tor  laboratory  animat,  human  brain, 
and  dedicated  PET  mammography  applications. 

II.  MATERIALS  AND  METHODS 

We  have  used  ribbons  of  wavelength-shilling  libers  to 
determine  the  positions  of  gamma-ray  interactions  wiihin  thin 
LSO  scintillator  crystals.  Within  a  wvaelencth-shifting  fiber, 
fluorescent  dopant  molecules  absorb  short-wavelength 
incident  photons  and  then  isotropically  emu  secondary  longer- 
wavelength  photons.  One  can  use  the  fluorescent  re-emission 


process  effectively  to  inject  light  into  a  light  pipe  through  the 
sides  of  the  pipe,  since  a  fraction  of  the  re-emitted  light  is 
trapped  and  piped  along  the  fiber.  In  this  way,  the  fiber 
becomes  a  light  sensor,  with  short-wavelength  light  as  input 
and  long-wavelength  light  as  output.  An  array  of  parallel 
fibers  (i.e.  a  fiber  ribbon)  may  then  be  used  to  measure 
incident  light  profiles  in  one  dimension;  two  perpendicular 
ribbons  (on  opposite  faces  of  the  thinsciniillator  crystal)  can 
provide  a  two-dimensional  incident  light  profile  measurement. 
Gamma-ray  interactions  in  dense  scintillators  result  in  a  very 
nearly  point-like  light  source,  because  of  the  short  range  of 
secondary  photocapture  electrons.  For  polished  high-index 
scintillator  crystals,  only  light  emerging  nearly  normal  to  the 
crystal  surface  escapes  total  internal  reflection,  resulting  in  an 
emission  spot  with  diameter  comparable  to  the  crystal 
thickness.  We  earlier  described  this  technique  as  applied  to 
BGO  crystals  read  out  with  red  wavelength-shifting  fibers  [1J. 

In  the  present  work  we  have  focussed  on  the  fast,  high-Z, 
high-brightness  scintillator  Lutetium  Oxyorthosilicate  or  LSO 
[2].  We  have  performed  most  of  our  initial  measurements  on 
two  samples  of  LSO,  each  2mm  x  10mm  x  10mm.  which  we 
obtained  from  Schlumberger-Doll  Research  [3].  These 
crystals  were  hand-polished  on  their  large  surfaces,  and  all 
surfaces  not  coupled  to  fibers  were  painted  black.  Fibers  were 
coupled  to  crystals  with  a  high-index  epoxy  (n=l.65),  EpoTek 
OG127  [4).  For  readout  fibers  we  used  l -mm  diameter, 
round,  doubly-clad  wavelength-shifting  fibers  from  Kuraray 
Corp  [5];  we  obtained  our  best  spectral  matching  with  their 
spectral  type  Y-l  1.  Doubly-clad  fibers  have  a  core  refractive 
index  of  1.60,  an  inner  cladding  index  of  1.49.  and  an  outer 
cladding  index  of  1.42;  this  yields  a  numerical  aperture  of  0.74 
and  a  trapping  efficiency  of  about  7%  per  fiber  end  for  the  re- 
emitted  light  (the  remaining  86%  of  the  re-emitted  light  is  not 
trapped  and  emerges  through  the  sides  of  the  fiber  ribbon). 

Our  initial  tests  have  been  carried  out  with  the  apparatus 
illustrated  schematically  in  Figure  I.  We  form  a  pencil  beam 
of  51 1-keV  gamma  rays  by  requiring  the  coincident  detection 
of  two  back-to-back  annihilation  gammas  produced  by  a 
nearly  point-like  Na-22  positron  source.  The  lirst  ot  the  two 
gammas  strikes  a  small  1mm  x  3mm  x  1 5mm  LSO 
“triggering"  crystal  through  its  1mm  x  3mm  entry  face,  while 
the  partner  ray  goes  through  a  1mm  lead  collimator  to  the 
“target"  2mm  x  1 0mm  x  1 0mm  crystal.  The  “target”  crystal  is 
coupled  to  a  ribbon  of  10  parallel  1mm  diameter  round 
wavelength-shifting  fibers  and  to  a  trigger  PMT.  Since  the 
light  collected  by  the  trigger  PMT  is  proportional  to  the 
energy  deposited  in  the  target  crystal,  we  are  able  to  select 
events  with  photocapture  of  a  pencil-beam  gamma  by  the 
target  crystal,  by  requiring  a  large  trigger  signal  in  coincidence 
with  a  signal  from  the  “triggering"  crystal.  The  “triggering" 
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LSO  crystal  and  its  PMT.  the  Na22  source  and  collimator  are 
all  mounted  on  a  calibrated  X-Y  adjustable  table,  while  the 
“target"  crystal  and  PMT  with  Tiber  assembly  are  fixed:  in  this 
way.  a  pencil  beam  can  be  precisely  directed  into  the  target 
setup,  and  may  be  swept  along  the  target.  AH  measurements 
are  made  with  unshaped  signals  directly  coupled  to  charge- 
sensitive  ADCs  (LeCrov  2249W);  the  timing  gate  width  was 
I50nscc  and  only  photopeak  events  were  used  tor  analysis. 


physics,  which  allows  for  simple  implementation  and 
modifications  of  the  detector  geometry  through  a  flexible  and 
powerful  user  interface.  We  generate  gamma-rays  from  our 
source  phantom  according  to  a  positron  range  distribution 

appropriate  to  18F.  and  with  acollinearity  between  the 
annihilation  gammas  sampled  from  a  Gaussian  distribution 
with  width  0.3  degrees.  The  detector  resolution  was  usually 
taken  to  match  the  (not-yet-optimized)  response  functions 
measured  with  our  test  apparatus,  and  matches  the  measured 
dependence  on  light  yield  (gamma  ray  energy). 

We  have  also  begun  work  on  making  images  with 
simulated  data  using  a  simple  2-D  filtered  back-projection 
reconstruction  algorithm,  for  a  number  of  test  phantoms.  In 
particular,  we  have  looked  at  the  benefits  of  depth-of- 
interaction  resolution  for  small-diameter  systems  such  as  for 
use  with  laboratory  animals.  For  some  of  these  studies  we 
have  simulated  a  1 /2-scale  version  of  a  standard  resolution 
phantom,  with  resulting  10cm  diameter  axially  symmetric 
phantom  viewed  by  a  simulated  20cm  diameter  detector  ring. 
To  contrast  with  our  detector,  we  have  modelled  a  device 
without  depth-of-interaction  sensitivity,  consisting  of  8mm 
thick  LSO  on  a  20cm  diameter  ring.  For  the  radially- 
segmented  fiber  readout  20cm  diameter  ring  we  have  used  a 
total  crystal  thickness  of  24mm.  which  would  result  in 
improve  coincident  efficiency  by  a  factor  of  3.2  relative  to  a 
8mm  thickness  ring. 

A  wide  range  of  system  design  options  are  open  with  a 
radially-segmented  fiber  readout  geometry,  and  we  have 
begun  our  exploration  of  these  options  with  a  baseline  module 
having  thickness  24mm  and  area  1 28mm  x  1 28mm.  Our  plan 
is  to  construct  and  test  at  least  one  such  module  within  the 
coming  year.  Assuming  2mm  thick  LSO  layers  and  1mm 
diameter  fibers,  this  requires  the  readout  of  13  x  128  =  1664 
fibers  per  module.  Multiplexing  is  clearly  necessary  in  order 
minimize  the  cost  for  readout  photomultipliers  and  digitizing 


Each  ot  10  fibers  coupled  to  the  target  crystal  is  read  out 
by  an  individual  Hamamatsu  R1635  PMT.  allowing  us  to 
record  the  light  profile  and  centroid  position  within  the  fiber 
ribbon  on  an  event-by-event  basis:  all  PMTs  were  calibrated 
and  outputs  convened  to  photoelectrons.  The  index  ot 
refraction  of  LSO  is  about  1.8.  so  we  used  the  high-index 
optical  epoxv  for  coupling  crystals  to  libers.  The  thin  sides 
were  blackened  with  an  opaque  mixture  of  high-index  epoxy 
and  xerographic  toner  (to  simulate  a  larger  crystal).  Fibers 
were  coupled  to  both  ol  the  large  LSO  surtaccs.  with  30cm 
long  fibers  having  one  polished  end  at  the  readout  PMTs  and 
their  other  ends  I  cm  from  the  crystal.  Light  trapped  at  the 
optical  interface  between. cladding  and  air  ("cladding  light  ) 
was  removed  by  an  opaque  fiber  holder  gnpping  each  fiber 
before  its  coupling  to  us  readout  PMT.  The  nber/PMT  optical 
intertacc  was  made  with  optical  grease. 

We  have  developed  a  complete  Monte  Carlo  simulation 
of  our  detector's  optics,  including  all  optical  interfaces  and  the 
correct  geometry  tor  round  multiclad  fibers.  We  have  also 
developed  a  detailed  radiation  transport  Monte  Carlo 
simulation  to  study  the  effects  of  gamma-rav  scattering  and 
attenuation  both  in  objects  being  imaged  and  within  the 
detector  itself.  The  latter  simulation  is  based  on  the  GEANT 
detector  simulation  code  commonly  used  in  high-energy 


electronics.  One  method  for  such  multiplexing,  employing 
separate  readout  of  two  fiber  ends,  is  illustrated  schematically 
in  Figure  2  below.  Other  schemes,  including  mirrored  fibers 
and  two  different  colors  of  waveshifting  fiber,  are  possible. 


One  Plane  ol  8  Parallel  Fiber  Ribbons 


Profile  across  eacn  ol  Index  tor  each  of 

8  fiber  nbbons  8  titoor  ribbons 

multiplexed  onto  1 6  anodes  f »wh  8  anodes 


Figure  2:  Multiplexed  readout  of  fiber  Ribbons  using  16- 
channel  multianode  photomultipliers. 
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We  have  been  testing  a  new  R5800-L16  multianode 
photomultiplier  from  Hamamatsu  Corporation,  which 
combines  16  input  pixels  each  I  mm  x  16mm  in  area  with  a 
very  narrow  pulses  (2ns  FWHM)  and  good  •ongle- 
photoelectron  resolution  |6).  By  reading  out  the  iwo  ends  ot 
muluclad  liber  ribbons  to  two  different  multianode  PMTs. 
one  can  read  out  as  many  as  256  liber  ends  per  PMT. 
Multiplexed  readout  of  fiber  ribbons  corresponding  to 
different  crystal  layers  is  also  possible  and  probablv  desirable. 
If  0.5mm-diameter  libers  are  used,  the  phxsical  limit  is  iOOO 
liber  ends  per  PMT.  Although  not  indicated  in  the  figure,  it 
may  be  desirable  to  break  up  each  128mm  x  1 28mm  crystal 
plane  into  64  16mm  x  16mm  optically  isolated  "cells"  to 
reduce  multiplexing  cross-talk. 

Electronics  channel  count  can  be  minimized  by  charge- 
division  readout  (two  signals  rather  than  16  irom  each 
photomultiplier)  as  is  typically  done  wuh  ^rosscd-wire 
position  sensitive  PMTs.  We  have  measured  .Ons  pulse 
widths  unto  a  50  ohm  load)  and  «ltnm  resolution  at  20 
photoelccirons  for  charge  division  readout  ot  our  R5300-L16 
PMT.  using  100  ohm  resistors  to  connect  the  anodes.  In 
principle  the  readout  of  a  16-laver  module  with  area  12S  mm 
x  123mm  using  0.5mm  libers  could  require  ijm  4  multianode 
photomultipliers  tat  a  current  cost  ot  SI500  cuchi.  with  8 
associated  fast  ADC  channels. 


in.  preliminary  results 

Our  measured  light  yield  tor  fiber  readout  ot  LSO  crystals 
is  in  general  agreement  with  our  expectations  The  light 
output'  for  LSChs  quoted  at  venous  levels  in  the  literature,  and 
is  a  strong  function  of  the  concentration  ot  Cerium  dopant:  the 
most  recent  work  cites  a  light  yield  ot  ).t mo  photons/MeV 
deposited  or  about  10.000  photons  per  annihilation  gamma 
photocapture  event  lor  LSO  with  a  Cerium  concentration  ol 
0.22Cr  |7|  Monte  Carlo  simulations  predict  that  more  than 
25?r  ot  emitted  photons  should  emerge  trom  a  given  crystal 
readout  tacc  into  the  corresponding  liber  ribbon,  with  photons 
incident  at  targe  angles  (relative  to  the  surtaec  normal)  totally 
internallv  reflected.  Direct  coupling  of  a  photomultiplier  with 
209c  quantum  efficiency  (when  averaged  over  the  LSO 
spectrum)  to  LSO  with  this  brightness  should  thus  give  about 
500  photoeiectrons:  this  is  consistent  with  our  measured  light 
vield  in  this  "direct  coupled"  geometry  The  next  link  in  the 
optical  chain  is  the  match  of  the  emission  spectrum  ot  the 
LSO  crvstal  to  the  absorption  spectrum  <<t  the  liber.  This 
tactor  is  about  60*3-.  based  on  convolution  ot  published  LSO 
emission  spectra  and  Kuraray  s  absorption  spectra  tor  Y-l  1. 
We  thus  expect  that  roughly  1500  scintillation  photons  are 
absorbed  and  rc-emnted  within  the  wavelength  shifting  libers 
in  each  photocapture  event. 

Muluclad  fibers  trap  14%  ot  their  rc-cmittcd  light,  or  7 9c 
per  fiber  end.  Our  measurements  comparing  light  yield  with 
muluclad  vs.  singly  clad  fibers,  where  the  latter  have  4 9c 
trapping  per  end',  confirms  this  trapping  traction  for  the 
muluclad  fibers  we  used.  This  implies  that  with  muluclad 
libers  about  1 00  photons  are  piped  toward  cacti  liber  end  each 
photocapture  event.  Finally,  some  photons  were  attenuated 
within  the  readout  libers  (such  self-absorption  lilters  out  the 


shortest-wavelength  photons,  and  typically  results  in  re- 
emission  of  an  unpiped  green  photon),  and  a  few  photons 
arriving  at  the  photomultiplier  were  reflected  at  the  optical 
interface  between  fiber  and  photomultiplier.  Our  estimates  of 
these  effects  were  that  we  had  about  80  photons  entering  our 
liber  readout  photomultipliers  each  photocapture  event,  where 
they  had  an  average  quantum  efficiency  of  about  10%  for 
being  convened  to  photoeiectrons  (Y-l  1  has  its  emission  peak 
at  495  nm).  This  accounts  for  the  roughly  8  photoeiectrons 
which  we  observed  per  photocapture  event  with  our  test 
apparatus.  This  light  yield,  although  small,  is  sufficient  for 
multiplexed  readout  with  good  efficiency.  Additional  light 
yield  may  be  achievable  by  using  fibers  with  a  better  spectral 
match  to  LSO  emission,  by  using  PMTs  with  better  quantum 
efficiency  for  green  light  (green-extended  R5800-L16  PMTs 
are  currently  under  development  at  Hamamatsu)  or  by 
obtaining  brighter  LSO. 

The  emission  light  profile  from  LSO  crystals  and  the 
measured  spatial  resolution  achieved  at  our  current  light  level 
is  shown  in  Figure  3  below.  The  resolution  horizontal  scale 
has  been  calibrated  by  moving  the  pencil  beam  in  our  test 
fixture  across  the  crystal  by  known  increments  and  measuring 
the  resulting  changes  in  the  measured  centroid  position. 


Figure  3:  Measured  response  of  LSO/Fiber  test  apparatus. 

Left  =  Predicted  and  measured  profile  of  light  collection. 

Right=  Measured  spatial  resolution  using  centroids. 

By  measuring  the  distribution  of  photoeiectrons  across 
the  fiber  ribbon  for  photocapture  events,  we  were  able  to 
reconstruct  the  event  position  with  an  accuracy  of  about 
0.7mm  RMS.  There  is  still  a  contribution  from  the  beam 
width  to  this  resolution,  for  which  we  have  not  corrected.  The 
accuracy  of  our  coordinate  measurement  varies  directly  with 
the  width  of  the  light  emission  profile,  so  that  the  quality  ot 
crystal  surface  finish  is  important.  Our  Monte  Carlo  optical 
simulation  predicted  a  profile  of  about  1 .0mm  RMS  for  a 
2mm  thick  ideally  polished  LSO  crystal,  while  the  measured 
profile  had  a  width  of  1.6mm  RMS.  Evidently,  our  crystal 
surface  preparation  and  fiber  coupling  could  be  improved. 
Spatial  resolution  improves  with  light  yield  like  the  square 
root  of  the  number  of  photoeiectrons  N,  and  ideally  the 
resolution  should  approach  the  profile  width  /  V  N. 
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We  have  carried  out  a  simulation  tor  a  geometry 
appropriate  to  Positron-Emission  Mammography,  which  we 
have  simplified  as  two  parallel  planes  separated  by  10cm  of 
water.  Figure  4  illustrates  the  line-spread  function  (LSF) 
predicted  for  a  perfect  detector  (due  to  just  positron  range  and 
acol linearity),  for  fiber  readout  trom  1mm  thick  LSO.  and  tor 
fiber  readout  from  2mm  thick  LSO:  the  latter  two  are 
extrapolated  from  our  non-optimized  test  results.  It  is 
apparent  that  range  and  acollineanty  etfects  are  \cry  small . 


Figure  4  Simulated  line-spread  Junction  and  intrinsic  limits. 

As  mentioned  earlier,  we  have  contrasted  the  simulated 
response  of  our  proposed  detector  with  that  ot  a  device  with 
comparable  spatial  resolution  but  without  depth -ot-interactiori 
sensitivity:  both  were  formed  into  20cm  diameter  rings. 


Figure  5:  Resolution  vs.  radial  position. 

TOP  =  No  DOI.  MID  =  2mm/  fiber.  BOT=  i  mm/  fiber 

Figure  5  shows  the  reconstructed  resolution  as  a  function 
ot  distance  from  the  center  of  the  field-ol-view  for  each 
device).  Again,  we  have  used  the  measured  but  not-yet- 
opumized  resolution  from  our  test  apparatus  as  input  to  this 
simulation:  in  addition,  we  have  indicated  the  resolution 
which  should  be  achievable  with  1mm  thick  LSO  layers  with 


our  fiber  readout.  While  our  device  has  very  nearly  uniform 
resolution  over  the  entire  field  of  view,  the  device  without 
depth-of-interaction  sensitivity  has  clear  resolution 
degradation  away  from  the  center  of  the  field  of  view. 

IV.  DISCUSSION 

Radially  segmented  readout  of  thin  LSO  crystals  with 
wavelength-shifting  fibers  has  the  following  advantages: 

•  Very  High  Spatial  Resolution 

•  Depth-of-Interaction  Measurement 

•  High  Rate  Capability 

•  Potential  for  Low  Cost 

•  Flexible  Readout  Options 

Unlike  most  designs  for  very  high-resolution  PET  detectors, 
this  design  scales  readily  to  large  systems  while  preserving 
high  efficiency  and  fine  resolution;  it  also  only  makes  use  of 
existing  photosensor  and  electronics  technology. 

Since  the  spatial  resolution  scales  linearly  with  the  crystal 
thickness,  extremely  fine  spatial  resolution  should  be 
achievable  with  small-diameter  rings.  Resolution  blurring  due 
to  acollinearity  (about  2.2mm  FWHM  per  meter  of  ring 
diameter)  can  be  made  to  dominate  system  resolution  for  large 
rings:  smaller  rings  can  achieve  submillimeter  resolution  with 
high  efficiency  across  a  large  Held  of  view  (if  short-range 
positrons  like  I8F  are  used)  by  using  very  thin  crystal  layers. 
Axial  resolution  is  naturally  comparable  to  tangential 
resolution,  and  to  depth-of-interaction  resolution.  Double-hits 
from  Compton-scatters  within  the  detector  ( which  can  degrade 
resolution  of  block  detectors)  are  readily  identifiable  with 
suitable  electronics,  since  they  will  almost  always  deposit 
energy  in  more  than  one  crystal  layer.  Use  may  be  made  of 
Compton  scattered  (rather  than  photocapturcd)  events  in  small 
rings  imaging  objects  with  little  scatter,  with  only  a  small 
degradation  in  spatial  resolution.  Finally,  detector  modules 
may  be  combined  in  a  "gapless"  geometry  by  radially 
overlapping  in  a  pinwheel  arrangement  as  shown  in  Figure  5: 
any  oversampling  which  which  results  in  the  overlap  regions 
may  be  corrected  using  depth-of-interaction  information. 


200  mm 


Brain  Ring 

Figure  5:  Conceptual  System  Geometries  with  Fiber  Readout 
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Depih-of-interacuon  measurement  is  critical  tor 
preserving  very  fine  spatial  resolution  throughout  a  large  field 
of  view,  without  compromising  on  detector  coincidence 
efficiency  (which  is  critical  for  the  statistics  needed  to  made  a 
high-resolution  image).  Depth-of-interaction  sensitivity 
permits  a  reduction  m  ring  diameter,  which  pays  several 
dividends:  lower  system  costs  tor  less  scintillator  and  fewer 
readout  channels,  smaller  effects  from  acollinearity.  and 
higher  efficiency  through  greater  solid  angle  coverage.  The 
last  comes  at  the  price  of  increased  scatter  acceptance,  but  this 
trade-off  is  already  being  chosen  in  the  selection  of  3-D  over 
2-D  PET.  Increasing  the  module  size  (both  in  axial  and 
tangential  dimension)  yields  a  more  cost-etfecuve  readout  by 
using  longer  readout  fibers,  which  increases  the  detector 
volume/PMT  area  rauo:  the  price  paid  is  lower  rate  capability. 

With  fully  multiplexed  readout  (including  in  depth)  of 
128mm  x  128mm  area  modules,  a  rate  capability  of 
>  I  Mhz/module  should  be  achievable.  The  45ns  decay  time  of 
LSO  permits  signal  integration  over  1 50ns.  with  digization  of 
charge-division  signals  to  the  required  7-bit  resolution 
possible  in  less  than  I  microsecond.  Time  correspondence  of 
8  PE  signals  from  45ns  scintillator  has  high  eiftcicncy  within 
a  I0-20ns  window.  The  coincidence  timing  tfor  singles 
rejection)  and  energy  resolution  achievable  with  an  all-fiber 
readout  are  more  problematic.  If  necessary,  fiber  readout  can 
be  supplemented  with  large  trigger  PMTs  which  would  collect 
the  re-emitted  but  non-trapped  light  which  emerges  through 
the  sides  of  the  wavelength-shifting  libers  and  to  which  both 
fibers  and  crystals  are  transparent.  Our  measurements  have 
shown  that  although  less  than  10%  ol  LSO  light  passing 
through  one  fiber  layer  is  absorbed  in  the  next,  the  layers 
which  are  opaque  to  LSO  light  are  very  transparent  to  fiber 
light.  Supplemental  readout  with  large  trigger  PMTs  would 
yield  hundreds  of  photoelectrons.  >utficient  lor  Ins 
coincidences  and  energy  resolution  limited  bv  etfects  other 
than  photosiatistics.  Large-area  photodiodes  might  be 
substituted  for  photomultipliers  to  preserve  compatability  with 
operation  in  a  magnetic  field,  if  necessary 

Fiber  readout  of  LSO  is  potentially  quite  cost-effective, 
presuming  that  the  bulk  price  of  LSO  continues  to  fall.  It  is 
currently  at  $50/cc  in  quantities  of  I00*s  ot  cc's.  and  is 
projected  by  the  crystal  grower  to  fall  to  SI5/cc  in  about  l 
year.  At  this  price  it  becomes  competitive  with  BGO,  and 
fiber  cutting  and  polishing  costs  become  signticant  (especially 
for  very  high-resolution  devices  with  small  crystals).  The 
crystal  cutting  and  polishing  required  for  the  thin  fiat  crystals 
used  with  fiber  readout  is  significantly  less  than  lor  the  long 
thin  crystals  in  other  detectors,  and  we  will  carefully  check 
how  our  performance  vanes  with  crystal  surface  treatment. 
The  cost  for  4  multianode  PMTs  to  read  out  fibers  from  a 
1 28mm  x  1 28mm  module  is  comparable  io  the  25  I"  PMTs 
which  would  be  required  to  read  out  the  same  area  in  a  block 
detector:  if  larger  modules  have  sufficient  rate  capability, 
fiber  readout  costs  fall  still  further.  The  cost  of  fibers  in 
quantity,  about  S  1/m  tor  0.5mm  diameter,  makes  the  fiber  cost 
small  in  comparison  to  crystals  and  PMTs.  Finally,  the 
electronics  channel  count  reduction  yields  savings,  reducing 
the  25  ADCs  for  the  above  block  detector  to  8  in  our  design. 


Although  our  initial  results  have  been  encouraging,  we 
have  been  somewhat  handicapped  by  the  availability  of  only 
small  samples  of  LSO.  This  situation  has  now  changed,  and 
we  anticipate  the  availability  of  hundreds  of  cc’s  of  LSO  by 
early  1996.  Our  plan  is  to  further  optimize  and  to  construct 
large  modules  with  this  new  material  as  soon  as  possible. 

V.  CONCLUSIONS 

We  have  demonstrated  that  wavelength-shifting  fibers 
can  be  used  to  read  out  gamma-ray  interaction  positions 
within  thin  sheets  of  LSO  crystals,  with  high  efficiency.  The 
spatial  resolution  achievable  (FWHM)  is  less  than  the 
thickness  of  the  crystal  layer  being  read  out.  Consequently, 
by  constructing  a  detector  from  a  thick  stack  of  many  thin 
layers,  one  can  combine  very  fine  spatial  resolution  with 
depth  of  interaction  sensitivity  and  high  efficiency.  With 
optical  multiplexing  at  the  crystal  layers  and  at  the  readout 
photomultipliers,  a  cost-effective  design  with  extremely  high 
spatial  resolution  is  achievable;  this  design  is  readily  scaled  to 
large-diameter  systems.  During  the  coming  year,  we  plan  to 
constract  one  or  more  128mm  x  128mm  x  20mm  thick 
modules  to  test  and  optimize  this  design. 
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Abstract ~  We  have  constructed  and  tested 
prototype  PET  detector  modules  incorporating 
wavelength-shifting  fiber  ribbons  alternating  with 
thin  CsI(Na)  scintillator  plates  in  multilayer  stacks. 
Modules  with  cross-sectional  area  of  11cm  x  11cm 
and  with  9mm  of  total  crystal  thickness  have  been 
tested.  A  reconstructed  spatial  resolution  of  3  mm 
FWHM  has  been  obtained  for  a  line  source,  with  an 
energy  resolution  of  <20%  and  a  time  resolution  o  f 
25ns  FWHM.  Energy  measurements  and  coincidence 
timing  are  provided  by  4  photomultipliers  in  an 
Anger  array  on  one  11cm  x  11cm  module  face. 
Multiplexed  fiber  readout  has  been  implemented  with 
multianode  photomultipliers  coupled  to  several  fiber 
ribbons.  The  methods  and  resolution  achieved  are 
readily  extended  to  larger  systems  in  a 
straightforward  and  cost-effective  manner.  Results 
of  tests  with  several  phantoms  are  presented. 

I.  Introduction 

Perpendicular  ribbons  of  wavelength-shifting  fibers  may 
be  used  to  determine  the  position  of  gamma-ray  interactions 
within  thin  scintillator  crystals,  with  the  obtainable  spatial 
resolution  comparable  to  the  crystal  thickness  [1-2].  By 
building  a  multilayer  stack  of  alternating  thin,  flat  crystal 
layers  and  perpendicular  fiber  ribbons,  we  have  measured 
gamma-ray  interaction  energies  and  positions  in  3  dimensions 
with  good  accuracy  and  efficiency.  We  have  constructed  and 
characterized  such  a  device,  and  have  imaged  3  test  phantoms. 

Depth-of-interaction  sensitivity  is  needed  to  preserve  the 
spatial  resolution  of  PET  detectors  when  accepting  events 
whose  lines-of-response  intersect  detector  elements  at 
significantly  oblique  angles.  The  parallax  which  leads  to 
errors  without  depth-of-interaction  sensitivity  may,  however, 
be  used  to  provide  information  on  object  depth-of-field  in  a 
limited  angle  geometry,  for  a  device  with  depth-of-interaction 
measurement  capability.  This  is  particularly  significant  in  the 
design  of  devices  for  Positron  Emission  Mammography 
(PEM)  where  two  parallel  plate  detector  elements  above  and 
below  a  breast  may  be  used  to  detect  coincident  gammas  [3-4]. 

For  a  device  whose  field  of  view  is  comparable  to  the 
detector  separation,  events  with  significant  parallax  may  be 
acquired  and  used  to  resolve  features  at  different  depths-of-field 
between  the  two  plates.  Other  groups  have  used  back- 
projection  imaging  in  this  geometry,  while  we  have  used  3-D 
limited-angle  tomographic  reconstruction  incorporating  an 
iterative  reconstruction  algorithm  based  on  Inverse  Monte 
Carlo  methods  and  Image  Space  Reconstruction. 


Detector  Schematic 


Reflector 

Figure  1 :  Schematic  Illustration  of  detector  geometry.  Three  thin  layers  of 
CsI(Na)  scintillator  alternate  with  perpendicular  fiber  ribbons.  Edge  fibers 
are  used  to  measure  depth-of-interaction,  while  Anger  PMTs  view  the 
laminated  fiber/crystal  block  through  a  light  mixer;  the  block  walls  are  lined 
with  white  reflector. 

A  schematic  illustration  of  our  device  geometry  is  given 
in  Figure  1.  For  flat  polished  crystals,  total  internal  reflection 
guides  scintillation  light  emitted  at  shallow  angles  with 
respect  to  the  crystal  surface  toward  the  crystal  edges,  while 
light  emitted  nearly  normal  to  the  surface  exits  the  crystal 
above  and  below  the  gamma  ray  interaction  point.  This 
exiting  light  position  may  be  sensed  locally  by  wavelength- 
shifting  (fluorescent)  optical  fibers,  which  absorb  the  primary 
scintillation  photons  and  isotropically  re-emit  secondary 
photons  at  longer  wavelengths. 

Perpendicular  fiber  ribbons  (X-fibers  and  Y-fibers)  on 
opposite  sides  of  a  thin,  polished  crystal  layer  may  thus  be 
used  to  measure  gamma-ray  interaction  positions  in  two 
dimensions;  further  fibers  at  the  edges  of  crystal  layers  (Z- 
fibers)  may  be  used  to  sense  scintillation  light  which  was 
totally  internally  reflected  within  a  layer,  providing 
information  on  depth-of-interaction  within  a  multilayer  stack. 
Alternatively,  one  may  read  out  fiber  ribbons  from  different 


depths  within  a  multilayer  stack  with  different  photosensors, 
but  this  is  costlier  than  the  edge-fiber  decoding  method. 

Coupling  scintillator  crystals  to  wavelength-shifting 
fibers  is  much  cheaper  than  direct  fiber  coupling  with 
transparent  fibers,  because  of  the  smaller  photoconverter  (e.g. 
multianode  photomultiplier)  area  required  to  read  out  fibers 
from  a  given  area  of  scintillator.  With  wavelength-shifting 
fiber  readout,  the  fiber  readout  area  increases  in  proportion  to 
the  scintillator  perimeter,  while  for  transparent  fiber  readout 
the  fiber  readout  area  increases  in  proportion  to  the  scintillator 
area.  Within  a  wavelength-shifting  fiber,  a  small  fraction  of 
the  secondary  photons  emitted  at  angles  near  that  of  the  fiber 
axis  are  piped  to  fiber  readout  photosensors,  while  most 
emerge  transversely  from  the  fibers  and  traverse  both 
scintillator  and  fibers  within  a  multilayer  stack  transparently. 
This  unpiped  light  can  be  collected  with  a  standard  Anger  array 
on  the  module  surface  to  provide  energy  and  trigger 
information  as  well  as  approximate  interaction  coordinates. 
These  approximate  coordinates  can  then  be  used  to  de¬ 
multiplex  the  readout  of  several  fiber  ribbons  which  may  then 
be  read  out  with  the  same  position-sensitive  multianode 
photomultiplier.  Multiplexed  fiber  readout  is  important  for 
limiting  both  system  component  costs  and  electronics  channel 
counts. 

The  above  method  decouples  the  event  energy,  timing,  and 
triggering  functions,  which  are  performed  by  the  Anger 
system,  from  the  precise  event  position  determination  in  3 
dimensions,  which  is  performed  by  the  fiber  readout  system. 
This  is  important  since  the  light  yield  from  the  wavelength- 
shifting  fibers,  although  easily  sufficient  for  precise  coordinate 
determination  and  high  efficiency,  is  insufficient  for  precise 
event  energy  and  timing  determination. 

The  performance  achievable  with  wavelength-shifting  fiber 
readout  of  scintillator  crystals  depends  critically  upon  the 
choice  of  scintillator  material,  the  quality  of  crystal  surface 
treatment  and  related  optical  coupling  parameters,  and  finally 
upon  the  device  geometry  and  in  particular  the  crystal  layer 
thickness.  In  the  results  presented  below  we  have  used  CsI(Na) 
crystals  after  earlier  work  with  Nal(Tl)  and  LSO,  because 
CsI(Na)  being  less  hygroscopic  than  NaI(Tl)  is  easier  to  work 
with  and  because  of  its  lower  cost  and  greater  availability  than 
LSO.  Our  choice  of  layer  thickness  was  such  as  to  provide  a 
reasonably  efficient  multilayer  stack  with  just  a  few  layers. 
Our  initial  measurements  and  device  characterization,  reported 
here,  are  with  a  device  with  just  3  crystal  layers;  this  resulted 
from  our  limiting  the  number  of  layers  and  the  amount  of 
multiplexing  in  the  readout  so  as  to  study  in  detail  the  effects 
of  multiplexing  in  depth. 

II.  Materials  and  Methods 
A.  Detector  Configuration 

Measurements  were  carried  out  with  two  detector  modules 
operating  in  coincidence.  Each  module  contained  3  layers  of 


polished  CsI(Na)  scintillator,  with  each  layer  3mm  thick  and 
having  a  114mm  x  114mm  lateral  extent.  These  crystal  layers 
were  polished  first  with  a  succession  of  sandpapers,  then  with 
Aluminum  Oxide  powders  in  mineral  oil-based  slurries. 
Crystal  layers  were  alternated  with  perpendicular  wavelength 
shifting  fiber  ribbons  (Bicron  BCF91A,  1mm  diameter 
multiclad),  with  each  ribbon  112mm  in  width  and 
approximately  30cm  in  length.  Each  crystal  layer  was  coated 
with  a  high-index  (n=1.65)  resin  to  provide  provide  an  index 
matching  layer  between  the  crystal  (n=1.82)  and  the  outer 
cladding  of  fiber  ribbon  (n=1.42).  Acrylic  frames  with  4mm 
thickness  were  used  to  assemble  a  multilayer  stack  of 
scintillator  crystals  alternating  with  fiber  ribbons,  and  to  hold 
depth-encoding  Z  fibers  in  position  at  two  crystal  edges. 

After  assembly,  each  multilayer  stack  was  cast  in  clear 
epoxy  (Emerson  and  Cummins  Stycast  1267)  and  coupled  to  a 
1”  thick  acrylic  light  mixer  with  the  same  epoxy.  The  two 
X  and  Y  fiber  ends  opposite  the  fiber  readout  ends  woe 
optically  coupled  to  plane  mirrors  (aluminized  mylar)  to 
increase  the  number  of  photons  reaching  the  readout  ends.  The 
base  of  each  module  was  optically  coupled  to  a  diffuse  white 
reflector,  as  was  the  outside  of  each  acrylic  light  mixer  to 
increase  the  light  yeild  reaching  the  trigger  PMTs. 

Each  module  was  coupled  with  optical  grease  (Dow 
Corning  Q2-3067)  to  an  Anger-type  array  of  4  Phillips 
XP3697/PA  PMTs,  each  with  a  60mm  x  60mm  envelope. 
Linear  xlO  preamplifiers  (Phillips  NE5205AN)  were  mounted 
on  each  photomultiplier’s  base  to  boost  gains.  Detector 
modules,  the  Anger  PMTs,  and  fiber  readout  PMTs  were 
mechanically  mounted  in  a  light-tight  support  structure  with 
signal,  high-voltage,  and  low-voltage  feedthroughs. 

Each  112mm  fiber  ribbon  was  split  into  8  bands  each 
containing  14  fibers.  The  16  bands  of  X-fibers  and  the  16 
bands  of  Y-fibers  were  grouped  into  8  sets  of  4-fold 
multiplexed  bundles,  with  a  pair  of  bands  within  a  layer 
(lateral  multiplexing)  bundled  with  the  corresponding  pair  from 
another  layer  (depth  multiplexing).  Each  bundled  set  of  4 
bands  was  held  in  a  custom  Delrin  fiber  vise,  polished,  and 
coupled  with  optical  grease  to  a  Hamamatsu  R5900-L16 
multianode  photomultiplier.  The  14mm  bundle  width  ran 
across  the  discrete  anodes  of  the  multianode  PMT,  each  of 
which  corresponds  to  a  1mm  x  16mm  photocathode  region. 
The  16  anodes  of  each  multianode  were  then  read  out  with  ID 
charge  division:  the  sum  of  the  charge  division  signals  A1+A2 
indicated  the  light  yield  from  each  bundle,  while  the 
normalized  coordinate  (A2-A1)/(A2+A1)  measured  the  average 
lateral  position  of  light  across  the  bundle  [5]. 

B.  Data  Acquisition  Electronics 

The  signals  from  each  Anger  tube  were  split,  with  one 
copy  used  for  triggering  and  the  other  delayed  for  ADC 
measurement.  The  4  trigger  signals  from  each  detector  module 
were  analog  summed,  shaped,  and  discriminated;  coincidences 


between  discriminated  analog  sums  from  the  two  modules 
provided  triggers  which  generated  ADC  gates  and  initiated 
readout.  Additional  copies  of  each  discriminated  analog  sum 
were  sent  to  TDCs  (LeCroy  2228A)  for  timing  masurements. 

Copies  of  the  Anger  signals,  as  well  as  all  charge  division 
signals  from  the  fiber  readout  PMTs,  were  delayed  by  200ns 
and  digitized  within  1  microsecond  gates  by  LeCroy  2249 W 
ADCs.  A  total  of  2  TDC  channels  and  47  ADC  channels  were 
read  out  through  CAMAC  for  each  event.  Data  was  acquired 
with  a  PC  through  a  custom  PCI/CAMAC  interface  and 
controller  obtained  from  the  Budker  Institute  for  Nuclear 
Physics  in  Novosibirsk,  Russia.  This  interface  was  able  to 
collect  triggers  at  a  rate  of  several  hundred  Hz,  which  was 
sufficient  for  the  phantom  studies  described  below.  We  are 
currently  developing  a  higher-rate  interface  for  clinical  use. 

C.  Event  Reconstruction 

All  fiber  readout  PMTs  were  calibrated  to  give  the  number 
of  photoelectrons  collected  from  each  fiber  bundle  on  each 
event,  which  was  greatly  facilitated  by  the  excellent  single 
photoelectron  pulse  height  distribution  of  the  R5900-L16 
PMTs.  For  each  event,  the  X-  and  Y-  bundles  with  the  most 
light  were  identified,  and  the  Anger  and  depth  information  used 
to  demultiplex  and  identify  which  of  4  bands  within  a  bundle 
produced  the  signal.  The  Anger  system  spatial  resolution  was 
sufficient  such  that  the  demultiplexing  was  quite 
unambiguous.  The  charge  division  information  from  the 
selected  X-  and  Y-  fiber  readout  PMTs  were  then  converted  to 
coordinates  across  each  fiber  ribbon  and  combined  with  the 
position  of  each  fiber  ribbon  within  the  apparatus  to  give 
global  X-,  Y-  and  Z-  coordinates. 

Analytic  single-pass  3D  reconstruction  is  impossible  for 
the  limited-angle  geometry  of  two  parallel  detector  plates  (as 
required  by  the  physical  contraints  of  PEM),  since  the 
geometry  does  not  satisfy  Orlov’s  criterion.  Nonetheless,  3D 
information  is  present  in  the  parallax  of  oblique  rays,  and  can 
be  exploited  through  iterative  reconstruction  algorithms.  We 
have  used  an  Inverse  Monte  Carlo  algorithm  [6]  loosely  based 
on  ML/EM  but  operating  in  image  space,  and  have  used  all 
triggered  events[7].  The  essence  of  this  algorithm  is  that 
collected  events  are  back-projected  across  a  3D  voxel  array,  and 
this  back-projected  data  is  then  iteratively  compared  with  back- 
projected  “pseudo-data”  generated  by  Monte  Carlo  simulation 
of  the  hypothesized  source  distribution  at  each  step  [8].  This 
algorithm  rapidly  accomplishes  most  of  its  convergence,  and  5 
iterations  were  used  in  this  work.  Reconstruction  times  for 
these  datasets  were  comparable  to  acquisition  times  when 
performed  on  a  simple  single-processor  workstation. 


surrounded  by  a  uniform  background  as  illustrated  in  Figure  2. 
This  was  accomplished  by  filling  a  17cm  diameter  by  6cm 
height  cylinder  with  uniform  Ge-68  activity  distributed 
throughout  an  epoxy  base;  an  18mm  diameter  hollow  tube 
traversing  the  cylinder  accomodates  inserts  with  activity 
concentrations  10X  the  background  and  volumes  of  1.0  and 
0.5cc  each.  The  total  activity  was  10  |iC. 
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Figure  2:  Hot-spot  phantom  with  two  inserts  surrounded  by  uniform 
background  activity.  The  activity  concentration  of  the  hot  spots  was  lOx  the 
activity  concentration  of  the  uniform  background. 


The  second  phantom  we  imaged  was  of  our  own 
construction,  and  is  illustrated  in  Figure  3  below.  It  consisted 
of  a  3.5mm  thick  acrylic  piece  perforated  by  the  hole  pattern 
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Figure  3:  Hot-spot  resolution  phantom.  Spots  with  diameters  of  4mm,  3mm, 
2.5mm,  2mm  and  l  mm  have  separations  of  3  times  their  diameters. 


D .  Phantoms 

The  first  phantom  we  imaged,  which  was  manufactured  by 
North  American  Scientific  and  provided  for  our  use  by  PEM 
Technologies,  Inc.,  consisted  of  two  hot-spot  inserts 


shown,  sandwiched  between  two  0.25mm  thick  layers  of 
uniform  activity.  The  3.5mm  thick  insert  was  placed  within  a 
4mm  height  cylinder  and  filled  with  F-18  in  an  aqueous 
solution.  The  effect  is  to  produce  hot-spots  with  an 
activity^ackground  ration  of  8:1;  the  initial  activity  of  this 


phantom  was  30|iCi  distributed  in  a  volume  of  1 .3  cm  cubed. 
Both  the  first  and  second  phantoms  were  oriented  with  their 
narrowest  dimension  normal  to  the  two  detector  planes,  and 
with  their  geometrical  center  mid-way  between  the  two  planes. 

The  third  phantom  which  we  imaged  consisted  of  two 
hypodermic  needles  (8cm  length  x  <lmm  ID)  containing  F-18 
in  an  aqueous  solution,  mounted  perpendicular  to  each  other 
and  lying  within  two  parallel  planes  separated  by  14mm; 
these  planes  were  then  oriented  parallel  to  the  two  detector 
planes,  with  one  of  the  two  planes  passing  through  the  point 
mid-way  between  the  two  detector  planes.  The  initial  activity 
of  this  phantom  was  also  about  30(iCi. 

III.  Experimental  Results 


Figure  4:  Collected  light  distribution  (in  photoelectrons)  at  trigger  level  for 

the  analog  sum  of  4  Anger  PMTs  within  a  detector  module.  The  photopeak 
is  clearly  visible  and  easily  selected  by  a  simple  threshold  trigger  circuit. 


A.  Detector  Performance 

Detector  system  performance  was  measured  in  response  to 
the  third  phantom;  separate  measurements  tested  performance 
of  the  Anger  energy/trigger/timing  subsystem  and  of  the  fiber 
coordinate  measurement  subsystem.  The  Anger  PMTs  within 
each  module  were  calibrated  with  an  LED  to  determine  the 
approximate  number  of  photoelectrons  per  ADC  bin,  taking 
the  event-to-event  fluctuations  in  response  to  the  LED  as 
entirely  due  to  photostatistics.  Figure  4  below  shows  the 
energy  resolution  (i.e.  the  ADC  sum  for  the  4  Anger  PMTs 
within  a  module,  after  calibration)  and  the  Anger  system  light 
yield  in  photoelectrons.  The  energy  resolution  obtained  at 
trigger  level  for  those  events  which  generated  a  trigger  was 
better  than  20%,  including  nonuniformities  in  detector 
response  both  laterally  (e.g.  less  light  was  collected  from  the 
detector  comers)  and  in  depth.  It  should  be  noted  that 
information  on  event  positions  both  laterally  and  in  depth  is 
available  on  an  event-by-event  basis,  so  that  the  obtainable 
offline  energy  resolution  is  significantly  better  than  that 
shown.  The  low-energy  tail  of  our  trigger  senstivity  reflects 
the  limitations  of  the  very  simple  shaping  circuit  used.  The 
timing  performance  for  two  detectors  in  coincidence  is  shown 
in  Figure  5.  demonstrating  25ns  FWHM  time  resolution. 

The  fiber  coordinate  measurement  achieved  nearly  3mm 
FWHM  spatial  resolution  in  response  to  a  line  source  within 
the  third  (needle)  phantom,  with  the  long  tails  resulting  from 
the  presence  of  the  second  line  source  outside  the  focal  plane  as 
shown  in  figure  6.  Measurements  with  a  point  source  showed 
better  spatial  resolution,  approaching  2mm  FWHM  when 
restricted  to  one  small  section  of  one  module  coincident  with 
another  small  section  on  the  opposing  module.  We  attribute 
this  to  variations  in  crystal  surface  finish  and  some  fiber 
readout  calibration  errors.  The  average  fiber  light  yield  for  the 
fiber  band  with  most  photoelectrons  in  a  photopeak  event  was 
10  photoelectrons,  as  shown  in  Figure  7. 
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Figure  5:  Time  difference  (in  nanoseconds)  between  signals  from 
discriminated  anode  sums  for  two  coincident  detector  modules. 
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Figure  6:  Measured  line  spread  function  in  response  to  a  line  source 
spanning  nearly  the  width  of  the  field  of  view. 


Figure  7:  Measured  light  yield  distribution  (in  photoelectrons),  for  one  end  of 
one  fiber  bundle  as  read  out  through  one  multianode  photomultiplier. 


B.  Phantom  Images 


The  image  reconstructed  from  data  collected  with  the  hot¬ 
spot  phantom  illustrated  schematically  in  Figure  2  is  shown  in 
Figure  8  below.  Both  the  lcc  and  0.5cc  hot-spots  are  clearly 
distinguished  from  the  uniform  background,  with  this  image 
reconstructed  from  about  100,000  photopeak-photopeak 
events.  The  detector  separation  was  >10cm  due  to  some 
mechanical  constraints,  but  the  depth-of-field  resolution 
achieved  in  nonetheless  quite  good.  The  geometric  acceptance 
for  coincidence  events  falls  nearly  linearly  with  distance  from 
the  device  edge,  so  that  edges  and  comers  of  the  image  are 
more  susceptible  to  statistical  fluctuations.  Nonetheless,  with 
just  100K  events,  a  single  threshold  is  sufficient  to  eliminate 
all  activity  outside  the  hot-spots  throughout  the  field  of  view, 
as  shown  in  the  3-D  rendering  at  the  bottom  half  of  Figure  8. 


The  image  reconstructed  from  data  collected  with  the  hot¬ 
spot  resolution  phantom  illustrated  schematically  in  Figure  3 
is  shown  in  Figure  9  below.  The  4mm,  3mm,  and  2.5mm 
hot-spots  are  completely  resolved,  as  are  some  are  of  the  2mm 
hot-spots.  The  total  number  of  events  is  again  about  100K. 


k 


Figure  9:  Reconstructed  image  of  the  hot-spot  resolution  phantom  (See 
Figure  3).  The  phantom’s  hot-spot/background  activity  ratio  was  8:1. 

The  image  reconstructed  from  data  collected  with  the  two- 
needle  phantom  is  shown  in  Figure  10  below.  The  spacing 
between  the  two  detector  planes  was  9cm,  and  the  two  needles 
are  easily  resolved  in  depth.  The  number  of  photopeak- 
photopeak  coincident  events  used  was  about  70K. 


V.  Conclusions 


Separation  of  functions  between  energy/trigger/timing 
(with  the  Anger  array)  and  interaction  coordinate  measurements 
(with  the  wavelength-shifting  fiber  readout)  allows  their 
separate  optimization.  High  spatial  resolution  results  from 
small  single  crystal  layer  thickness,  while  high  efficiency 
results  from  a  large  total  crystal  thickness.  Separate 
optimization  thus  requires  the  capability  to  preserve  Anger 
system  performance  in  a  device  with  a  large  number  of  layers. 
The  initial  studies  reported  here  indicate  the  capability  of  3- 
layer  modules,  and  we  are  continuing  this  work  on  modules 
containing  more  crystal  layers.  Preliminary  measurements 
have  shown  good  energy  resolution  at  the  trigger  level  for  7- 
layer  modules,  when  adequate  care  is  taken  in  crystal  surface 
treatment.  The  quality  of  the  crystal  surface  polish  is  critical 
both  to  the  fiber  coordinate  measurement  and  to  the  Anger 
energy/trigger  function.  For  the  fiber  coordinate  measurement, 
it  determines  the  degree  to  which  photons  spreading  and  totally 
internally  reflecting  away  from  the  interaction  point  are  able  to 
exit  the  crystal  and  enter  fibers  before  reaching  the  crystal  edge. 
For  the  energy  measurement,  it  determines  the  overall 
transparency  of  a  multilayer  crystal/fiber  stack  to  wavelength- 
shifted  light  which  is  not  piped  down  fibers.  Although 
CsI(Na)  is  soft,  it  is  somewhat  hygroscopic  and  difficult  to 
polish  well,  which  we  believe  contributes  significantly  to  the 
limits  on  our  device’s  spatial  resolution. 

Other  scintillator  material  properties  which  affect  the 
detector  performance  are  its  density  and  its  photofraction.  A 
Compton  scatter  interaction  where  the  secondary  gamma  is 
convened  elsewhere  within  the  detector  will  produce  an  Anger 
signal  which  lies  within  the  photopeak.  Right-angle 
Compton  scatters  which  convert  in  the  same  crystal  plate  are 
indistinguishable  from  photocaptures,  but  result  in  distributed 
light  production  within  the  scintillator  crystal  and  thus 
positioning  errors.  Positioning  errors  due  to  detector 
Compton  scattering  are  lessened  when  using  a  scintillator  with 
a  higher  photofraction,  and  a  denser  scintillator  can  be  used  to 
decrease  the  number  of  layers  and/or  to  decrease  layer  thickness 
while  maintaining  detector  efficiency.  All  these  reasons,  in 
addition  to  its  greater  speed  relative  to  CsI(Na),  make  LSO  an 
excellent  material  for  very  high-performance  wavelength- 
shifting  fiber  readout  modules.  We  have  recently  acquired 
more  significant  amounts  of  this  material  and  are  proceeding  to 
develop  new  modules  containing  LSO  in  the  near  future. 

One  major  advantage  which  CsI(Na)  modules  maintain 
relative  to  LSO  modules  is  the  much  lower  cost  and  greater 
availability  of  CsI(Na).  For  that  reason,  we  are  also  planning 
to  construct  devices  with  significantly  larger  fields  of  view 
using  this  relatively  inexpensive  scintillator  material.  The 
cost  of  wavelength-shifting  fiber  readout  scales  proportional  to 
the  perimeter  of  the  device,  so  that  larger  areas  can  be  covered 
quite  cost-effectively.  Rate  capability  issues  are  significant  for 
larger  devices,  but  the  great  gain  in  efficiency  relative  to  thin 
crystals  (such  as  in  commercial  dual-plane  coincidence 
imagers)  shows  potential  for  significant  performance  gains. 


We  have  constructed  and  tested  two  multilayer  CsI(Na) 
PET  detector  modules  with  wavelength-shifting  fiber  readout, 
and  have  used  them  to  image  several  phantoms.  Wavelength- 
shifting  fiber  readout  has  intrinsic  depth-of-interaction 
sensitivity,  which  is  important  when  scintillators  are  brought 
into  close  proximity  to  an  object  being  imaged.  With  depth- 
of-interaction  sensitivity,  the  parallax  information  in  oblique 
incident  rays  can  be  used  to  resolve  objects  in  the  depth-of- 
field,  as  we  have  demonstrated.  In  a  two-plate  detector 
geometry  this  requires  iterative  reconstruction  techniques  and 
limited-angle  tomography;  and  we  have  shown  that  this  is 
feasible  for  a  device  applicable  to  Positron  Emission 
Mammography  (PEM).  Finally,  the  materials  used  to 
construct  our  test  device  were  relatively  inexpensive,  and  the 
device  design  scales  to  larger  sensitive  areas  cost-effectively. 
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Abstract 

We  have  developed  a  3D  PET  reconstruction  algorithm 
based  upon  Inverse  Monte  Carlo  analysis  and  have  tested  it 
on  data  acquired  by  high-resolution  3D  PET  detectors.  This 
algorithm  back-projects  3D  PET  data  across  a  3-dimensional 
array  of  voxels,  then  during  each  iteration  generates  Monte 
Carlo  “pseudo-data”  which  is  similarly  back-projected 
across  a  3-dimensional  voxel  array.  Comparison  between 
the  back-projection  of  the  data  and  the  back-projection  of 
the  pseudo-data  is  used  to  update  the  source  distribution 
hypothesis.  Inverse  Monte  Carlo  methods  provide  a  natural 
framework  for  incorporation  of  detector  response  functions, 
and  potentially  for  the  incorporation  of  effects  due  to  scatter 
and  other  sources  of  systematic  error.  The  algorithm  is  fairly 
fast,  capable  of  reconstructing  a  128  x  128  x  128  Hoffman  brain 
phantom  image  in  a  few  hours  on  a  modern  single-processor 
workstation.  We  have  achieved  nearly  linear  speed-up  by 
implementation  of  the  algorithm  on  a  computer  with  8  parallel 
processors. 

I.  Introduction 

The  recent  trend  in  operating  PET  detectors  in  three- 
dimensional  mode,  to  achieve  higher  sensitivity,  has  aroused 
interest  in  3D  image  reconstruction  algorithms.  The 
requirements  for  such  algorithms  is  that  they  have  to  be 
quantitatively  accurate,  as  well  as,  have  fast  reconstruction 
times.  The  3DRP  algorithm  proposed  by  Kinahan  and 
Rogers  [1]  is  most  commonly  implemented  on  commercial 
PET  scanners.  This  algorithm  is  robust  and  has  good  image 
contrast  but  suffers  from  poor  noise  properties  and  treatment 
of  statistical  data.  The  ML-EM  algorithm  [2,  3,  4],  widely 
applied  in  the  2D  tomography  is  computationally  cumbersome 
in  3D.  The  ML-EM  algorithm  is  an  iterative  algorithm,  which 
requires  calculation  of  the  system  response  matrix.This  is  the 
matrix  of  probabilities  that  a  photon  emitted  from  a  certain 
point  will  be  detected  in  a  particular  detector  line  of  response. 
This  probability  matrix  is  usually  too  huge  to  store  in  memory 
and  has  to  be  either  stored  in  on  disk  or  calculated  on  the  fly. 
Several  groups  have  used  symmetries  and  sparseness  of  this 
matrix  to  be  able  to  store  it  in  memory  [5].  Hudson  and  Larkin 
proposed  OS-EM  algorithm  [6]  which  uses  symmetries  in 
the  data  to  speed  up  the  iterations  but  still  quite  slow  in  3D. 
Gaining  popularity  for  3-D  reconstruction  are  the  approximate 
rebinning  algorithm  like  the  FORE  algorithm  which  rebins  data 
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into  set  of  2D  slices.  The  rebinned  slices  can  be  reconstructed 
using  filtered  backprojection  or  ML-EM,  OS-EM  or  WLS 
methods  [7]  .These  algorithms  are  computationally  fast  but 
the  accuracy'  of  the  reconstructed  images  is  limited  by  the 
approximations  implicit  in  the  line  integral  model. 

Our  interest  in  image  reconstruction  arises  from 
development  of  a  unique  detector  made  wavelength  shifting 
fibers  [8].  This  has  2  large  planar  detectors  facing  each  other 
the  position  coordinates  for  the  detected  events  are  readout 
by  1mm  fibers  and  our  detector  has  a  spatial  resolution  of 
2-3mm.  Design  considerations,  especially  for  mammogarphy, 
prohibits  construction  of  a  detector  geometry  with  complete 
angular  views.  This  incomplete  set  of  projections  excludes 
the  possibility  of  using  FBP  as  they  do  not  satisfy  Orlov’s 
criteria  [9].  ML-EM  algorithms  which  is  not  restricted  to 
device  geometry  is  a  possibility  but  the  computational  burden 
associated  with  explicitly  calculating  the  system  response 
function  makes  it  impartical.  Thus  we  felt  a  need  for  a  new 
algorithm  which  would  work  for  a  limited  angle  geometry 
detector  which  could  generate  fast  3D  images. 

Our  algorithm  is  based  on  Monte  Carlo  method  and  it  solves 
the  Inverse  problem  of  solving  for  the  source  distribution  given 
the  data  from  the  source  distribution.  Monte  Carlo  method 
for  solving  for  source  distribution  has  been  proposed  more 
than  10  years  ago  for  SPECT  imaging  [10].  Monte-Carlo 
methods  , however,  have  been  used  in  PET  to  study  the  effects 
of  attenuation  and  scatter. 

This  algorithm  is,  however,  not  restricted  to  a  dual  plane 
geometry  but  can  be  applied  to  a  wide  range  of  detector 
geometries.  We  will  present  its  implementation  on  a  two  high 
resolution  detectors  with  very  different  detector  geometries. 

II.  Theory 

The  Inverse  Monte  Carlo  algorithm  is  an  iterative  algorithm 
which  updates  a  hypothesized  source  distribution.  It  does  this 
by  Monte  Carlo  sampling  from  the  source  distribution  and 
comparing  it  with  the  acquired  data.  The  algorithm  digitizes 
the  image  space  into  voxels  and  starts  with  first  estimate  for  the 
source  distribution.Monte  Carlo  events  are  generated  from  the 
first  estimate  according  to  the  activity  of  the  source  distribution 
in  the  voxel.  In  other  words,  photon  pairs  are  generated  from 
a  from  a  random  point  within  the  voxel  and  given  a  random 
direction.  This  photon  pair  is  forward  projected  to  find  if 
they  hit  the  detector.  If  they  hit  the  detector  then  the  event  is 
backprojected  Backprojection  is  done  similar  to  the  method 
described  by  Barresi  [11]. 

This  algorithm  is  similar  to  the  image  space  reconstruction 
algorithm  (ISRA)  proposed  more  than  10  years  ago  by  Daube- 


Witherspoon  and  Muehllehner  [12]  in  that  the  measured  data  is 
stored  in  the  backprojected  array.  Also  comparison  of  the  data 
and  simulated  data  is  done  in  backprojection  space. 

The  difference  between  the  backprojection  of  the  data  and 
the  backprojection  of  the  simulated  data  is  used  to  update 
the  image  hypothesis.  In  the  following  iterations  the  number 
of  events  from  a  voxel  generated  is  equal  to  the  difference 
between  the  two  backprojections  at  that  voxel.  The  events 
that  are  subsequently  generated  are  similarly  projected  to  find 
where  they  hit  the  detector  and  then  backprojected. However,  if 
the  simulated  backprojection  at  the  source  voxel  is  less  than  the 
measured  backprojection  at  that  voxel  then  the  backprojection 
is  added.  In  the  case  that  the  simulated  backprojection  is 
more  than  the  measured  backprojection  the  backprojection  is 
subtracted. 

For  the  choice  of  first  estimate  for  source  distribution,  we 
can  start  with  a  uniform  source  or  image  formed  from  some 
earlier  iteration.  Also  we  can  also  take  the  backprojection  of 
the  data  as  the  first  estimate  and  it  has  been  seen  that  this  is  a 
very  good  starting  point. 

We  have  worked  with  a  stopping  criteria  that  the  number 
of  events  generated  is  about  10  times  the  number  of  events 
detected.  We  have  assumed  this  without  proof  but  we  have 
discovered  empirically  that  the  image  quality  does  not  improve 
significantly  for  more  events  generated. 


III.  Implementation  of  algorithm  for  2 

PLANE  GEOMETRY 


A.  Software  phantom  studies  for  BU  PET  geometry 


The  algorithm  was  tested  on  simulated  phantom  setup  with 
the  detector  model  based  on  the  detector  prototype  that  has 
been  fabricated  in  our  lab  [8].  The  detector  geometry  consisted 
of  2  parallel  plates  each  1 1cm  xl  1cm  ,  separated  by  6  cm  as 
illustrated  in  Fig  1.  This  detector  resolution  was  taken  to  be  2 
mm. 
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Figure  1:  Schematic  illustration  of  the  Boston  University  PET 
detector.Therc  are  only  2  planar  detectors  of  size  1 1  cm  x  1 1  cm  facing 
each  other.  The  distance  of  separtion  can  be  varied. 


We  simulated  20  million  events  with  an  array  of  64x64x32 


voxels  each  1mm  x  1mm  x  1mm  with  Poison  fluctuations.  We 
include  3  hot  spots  with  8  times  background  intensity;  these 
cylinders  are  10mm, 8mm  and  6mm  in  diameter  and  lengths 
1.5  times  their  diameter.  The  cylinders  will  be  reffered  to  as 
cylinders  A,B  and  C.  This  software  phantom  was  chosen  to 
match  that  used  in  a  previously  published  work  for  Positron 
Emission  Mammogarphy  [13].  Cylinders  A  and  B  lie  in  the 
same  Z  plane  and  can  be  seen  in  the  same  slice  as  shown  in  fig 
2(1,2). 


Figure  2:  (1)  3D  rendering  of  the  simulated  phantom  with  three 
cylinders- A,B,C.  (2)Slice  through  one  plane  showing  relative  intensity 
of  cylinders  A  and  B  versus  the  background  8:1.  Note  cylinder  C  is  on 
another  plane  and  cannot  be  seen.(3)3D  rendering  of  the  reconstructed 
image. (4)Slice  through  one  plane  of  reconstructed  image  showing 
relative  intensity  of  cylinders  A  and  B  versus  background. 


The  images  shown  in  Fig  2  (3,4)  was  reconstructed  in  25 
iterations.  The  convergence  of  the  algorithm  is  illustrated  in  fig 
3.  This  image  was  reconstructed  in  less  than  10  minutes  when 
implemented  in  a  SGI/Cray  Origin  2  Series  computer  with  a  192 
MIPS  R 10000  processor  running  at  195  MHz. 

5.  Phantom  studies  with  actual  data  from  BU  PET 

A  number  of  objects  were  imaged  and  reconstructed  using 
our  3  layer  CsI(Na)  prototype  detector.  This  work  has  been 
published  in  an  earlier  work  regarding  work  regarding  the 
performance  of  the  detector  and  so  will  not  be  reproduced 
here  [8]. 

IV.  Data  from  ring  detector  geometry 

This  algorithm  was  tested  on  data  taken  by  the  HEAD 
PENN-PET  detector  of  the  Hoffman  brain  phantom.  The 
HEAD  PENN-PET  scanner  is  a  volume-imaging  PET  scanner 
without  inter-plane  septa  with  high  spatial  resolution  and 
sensitivity  [14].  The  spatial  resolution  is  3.5  mm  in  both 
transverse  and  axial  directions  in  the  center  of  the  field  of 
view  (FOV).  The  FOV  is  25.6  cm  in  both  transverse  and  axial 
directions. 


A.  Format  of  the  data 


Figure  3:  The  Chi-Sq  convergence  rate  versus  number  of  iterations 


The  data  is  stored  into  four-dimensional  projection  matrices 
(x’,y\phi, theta)  that  are  128x128x96x15 .  The  spatial  sampling 
is  2mm.  The  data  had  been  normalized  for  efficiency  t  detector 
non-uniformities  and  attenuation  corrected. 


Symmetries  for  Cylindrical  Detector 


2  fold  symmeuy  in  2  direction 

8  fold  syraraeuy  in  x~y  plane 

Hence  accumulated  projections 
can  be  backpcojcctcd  as  a  set  of  16. 


Z 
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Figure  5:  Illustration  of  the  16  symmetric  lines  of  responses  for 
cylindrical  detector  geometry  and  rectangular  reconstruction  volume 


Figure  4:  Schematic  illustration  of  geometry  of  HEAD  PENN-PET 
detector  The  1 5  theta  projections  can  be  paired  into  8  sets  of  projection 
views  symmetric  about  the  XY  plane. 


5.  Symmetries 

The  four  dimensional  projection  matrix  which  is 
approximately  46  MBytes  is  too  huge  to  fit  into  memory  for 
most  computers.  Implementation  of  the  algorithm  for  this 
data  set  was  done  by  using  the  symmetries  presented  by  the 
cylindrical  geometry  of  the  detector.  The  projection  data 
was  divided  into  8  sets  of  projections  such  that  each  set  of 
projection  would  lie  within  2  symmetrical  projection  views 
about  the  XY  plane  as  illustrated  in  Fig  4.  The  algorithm 
was  modified  such  that  it  would  loop  through  the  8  projection 
views.  For  each  projection  view  the  events  generated  from 
the  source  distribution  were  confined  to  a  limited  theta  angle 
range  such  that  they  fell  within  that  projection  view.  The 
simulated  events  were  forward  projected  and  accumulated  in  a 
128  x  128  x  96  x2  projection  matrix.  For  each  projection  view 
backprojection  was  done  once  all  the  events  were  generated 
and  forward  projected  and  accumulated.  By  accumulating  the 
projection  and  backprojecting  simultaneously  we  could  use  the 
16  symmetries  in  the  lines  of  responses  Fig  5  to  give  a  linear 
speed  up. 

C.  Parallel  Implementation 

Since  the  projection  data  was  be  split  into  8  projection  views 
this  algorithm  could  easily  parallelized  by  8  processors  such 
that  one  processor  reconstructs  1  projection  view.  This  was 
implemented  on  a  a  SGI/Cray  Origin2  Series  parallel  computer 
system  with  multiple  processors.  Each  processor  was  a  192 
MIPS  R 10000  running  at  195MHz.  The  algorithm  was  written 


in  Fortran  and  the  parralelization  was  done  by  using  Message 
Passing  Interface  (MPI). 

D.  Results 

Data  the  Hoffman  Brain  Phantom  were  reconstructed  by  this 
algorithm  .  The  total  number  of  attenuation  corrected  events 
in  this  data  were  about  800  million.  Using  4  mm  size  voxels 
the  image  were  reconstructed  in  25  iterations  within  15  minutes 
while  simulating  3  billion  events, Fig  6.  Using  the  same  number 
of  iterations  and  simulated  events  but  using  2mm  voxel  size  the 
reconstruction  time  was  100  minutes.Fig7. 


Figure  6:  Hoffman  Brain  Phantom  using  4  mm  voxel  size;  using  25 
iterations  and  simulating  3  billion  events 


Figure  7:  Hoffman  Brain  Phantom  using  2  mm  voxel  size  ;using  25 
iterations  and  simulation  3  billion  events 


V.  Discussion 

The  images  from  the  2mm  voxel  size  had  higher  contrast 
but  were  noisier  that  produced  by  using  4  mm  voxel  size.  The 
image  quality  does  not  look  remarkedly  different  but  this  is 


due  to  the  effects  of  detector  resolution  which  is  3.5  mm.  To 
deconvolve  the  effects  of  detector  we  will  reconstruct  images 
from  the  digital  Hoffman  brain  phantom.  The  effects  of  the 
number  events  to  be  simulated  as  a  function  of  voxel  size  will 
be  studied  using  this  phantom  [15]. 

VI.  Conclusion 

We  have  explored  the  possibility  using  Inverse  Monte  Carlo 
methods  for  image  reconstruction  and  applied  it  to  data  from 
2  high  resolution  detectors  with  very  different  geometries. 
The  algorithm  is  shown  to  give  qualitatively  accurate  images 
of  complex  phantoms.  Precise  quantitative  analysis  of  the 
algorithm  will  be  done  next  using  the  EVAL3DPET  software 
developed  by  the  Medical  Image  Processing  group  at  the 
University  of  Pennsylvania  [16,  17]. 
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Abstract 

We  have  been  developing  3D  reconstruction  procedures  for 
detectors  with  a  fixed  dual-plane  acquisition  geometry,  and 
in  particular  for  Positron  Emission  Mammography  (PEM) 
detectors.  Because  this  limited-angle  acquisition  geometry 
does  not  satisfy  Orlov’s  criterion,  Fourier-based  analytic 
reconstruction  algorithms  may  not  be  used.  Nonetheless, 
depth-of-field  (tomographic)  information  is  present  in  the 
parallax  of  oblique  lines-of-response,  and  may  be  recovered 
using  3D  iterative  reconstruction  techniques.  The  mapping 
of  a  3D  source  volume  onto  a  4D  line-of-response  array 
requires  a  7D  system  response  matrix,  which  is  at  the  same 
time  inordinately  large  and  sparsely  populated.  We  have  taken 
advantage  of  translational  symmetries  for  detector  acquisition 
planes  parallel  to  planes  of  voxels,  and  have  stored  only 
non-zero  system  response  matrix  elements  as  a  list.  As  a 
result,  we  are  able  to  store  pre-calculated  system  response 
matrix  elements  for  a  36  x  36  x  40  voxel  array  sandwiched 
between  two  parallel  18x18  arrays  of  detector  crystals  in  less 
than  4  Megabytes  of  memory.  Because  these  matrix  elements 
easily  fit  in  computer  memory,  a  full  EM  iteration  requires  less 
than  10  CPU  minutes  on  a  single-processor  workstation.  We 
are  currently  examining  EM,  OSEM  and  ISRA  3D  iterative 
reconstruction  algorithms  using  this  compressed  system  matrix. 

In  addition  to  direct  application  of  the  system  matrix 
during  forward-  and  back-projection  steps  in  each  of  the  above 
algorithms,  it  is  possible  to  approximate  these  steps  using 
Monte  Carlo  sampling  techniques.  These  ’’Inverse  Monte 
Carlo**  methods  use  an  implicit  form  of  the  system  response 
matrix  contained  within  the  Monte  Carlo  simulation  of  the 
source  and  detector  systems.  In  the  context  of  a  dual  plane 
acquisition  geometry  we  are  therefore  able  to  compare  the 
speed  and  accuracy  of  Inverse  Monte  Carlo  and  direct  system 
matrix  application  techniques,  in  particular  as  a  function  of  the 
number  of  detected  events.  Results  of  quantitative  performance 
measures  using  simple  software  phantoms  are  presented,  as 
well  as  results  using  data  from  real  detectors  and  phantoms. 

I.  Introduction 

Iterative  reconstruction  methods  are  widely  used  in 
Emission  Computed  Tomography  (ECT)  to  more  accurately 
incorporate  source  and  detector  system  response  information, 
and  to  more  accurately  treat  event  statistics,  than  is  possible 

*This  research  was  supported  by  US  Army  Medical  Research  and 
Material  Command  Breast  Cancer  Research  Program  Grant  DAMD17- 
96-1-6192  and  in  part  by  the  Commonwealth  of  Massachusetts 
Department  of  Public  Health  Breast  Cancer  Research  Program  Grant 
SC- DPH- 3408-799 D049. 


with  less  computationally  burdensome  analytic  reconstruction 
methods  [1].  Iterative  techniques  have  been  much  more 
widely  used  in  2D  rather  than  3D  acquisition  geometries  [2], 
because  the  resulting  system  response  matrix  is  4-dimensional 
for  2D  acquisition  (a  2D  array  of  lines-of-response  for  each 
source  element)  and  7-dimensional  for  3D  acquisition  (a  4D 
array  of  lines-of-response  for  each  source  element).  Direct 
storage  of  a  7D  array  with  any  reasonable  number  of  voxels 
and  lines-of-response  is  usually  impractical,  even  after  array 
compression  to  save  only  non-zero  elements  and  application 
of  system  symmetry  operations  [3],  Recently,  methods  have 
been  developed  to  rebin  3D  data  as  a  set  of  nearly  equivalent 
2D  projections,  but  these  operations  inevitably  trade  off 
some  spatial  and  statistical  accuracy  for  the  considerable 
computational  advantages  of  parallel  2D  rather  than  fully  3D 
reconstruction.^] 

We  have  been  interested  in  the  application  of  Monte  Carlo 
sampling  techniques  to  the  forward-  and  back-projection 
operations  required  in  3D  iterative  reconstruction  algorithms, 
where  the  information  which  might  be  stored  explicitly  in  the 
system  response  matrix  is  embedded  implicitly  in  the  Monte 
Carlo  system  simulation  model.  This  allows  for  fully  3D 
reconstruction  with  no  re-binning  and  accurate  modelling 
of  the  statistical  properties  of  the  data.  Such  methods  have 
the  further  attractive  feature  of  potentially  incorporating  a 
detailed  model  of  the  source  (e.g.  attenuation,  scatter)  and 
detector  system  (e.g.  resolution,  depth-of-interaction,  detailed 
geometry)  properties.  The  computational  burden  associated 
with  such  modelling  is  less  than  it  might  at  first  appear  (when 
efficiently  coded),  and  it  becomes  less  burdensome  as  more 
computer  power  and/or  parallel  processing  capability  becomes 
available. 

Our  group’s  interest  in  3D  iterative  reconstruction 
techniques  grew  out  of  our  need  for  tomographic  image 
reconstruction  capability  in  a  limited-angle  data  acquisition 
geometry,  in  particular  for  application  to  a  Positron  Emission 
Mammography  (PEM)  detector  [5].  This  device  has  a  fixed 
dual-plane  acquisition  geometry,  acquiring  depth-of-field 
information  along  the  axis  between  the  two  detector  planes  in 
the  form  of  parallax  between  views  along  different  oblique 
lines-of-response.  Analytic  reconstruction  techniques  are  ruled 
out  because  this  limited-angle  acquisition  geometry  violates 
Orlov’s  criterion  [6].  We  were,  however,  able  to  develop  a  3D 
reconstruction  procedure  using  a  form  of  Inverse  Monte  Carlo 
analysis,  and  in  particular  through  a  Monte  Carlo  sampling 
implementation  of  the  ISRA  (Image  Space  Reconstruction 
Algorithm  [7])  method.  Preliminary  results  from  this  Inverse 
Monte  Carlo  approach  are  currently  under  review  [8]. 


Recently,  we  have  been  developing  a  very  compact 
dual-plane  PEM  detector  for  application  to  image-guided 
biopsy,  consisting  of  two  parallel  planes  each  with  an  18  x 
18  array  of  GSO  crystals,  and  with  each  crystal  3mm  x  3mm 
x  15mm.  In  the  process  of  developing  image  reconstruction 
methods  for  this  new  device,  it  became  apparent  that  the 
dual-plane  acquisition  geometry  has  sufficient  symmetry  to 
permit  storage  of  a  pre-calculated  system  response  matrix  in  a 
strikingly  small  amount  of  computer  memory.  Consequently, 
this  has  allowed  us  to  evaluate  the  merits  and  limitations  of 
Monte  Carlo  sampling  calculations  as  opposed  to  direct  system 
matrix  element  calculations,  by  direct  comparison  of  both 
techniques  as  applied  to  simulations  and  measurements  with 
our  new  PEM  device.  In  what  follows  we  will  present  some 
first  results  from  these  comparative  studies. 

II.  Theory 

We  use  the  notation  of  Titterington  [9],  as  have  others. 
For  J  source  pixels,  the  jth  element  has  emission  density 
A j.  The  measured  data  n  *  is  the  number  of  coincidences 
recorded  in  the  ith  of  /  pairs  of  coincident  detector  elements. 
The  conditional  probability  that  an  event  emitted  from  source 
element  j  is  assigned  to  projection  i  is  given  by  (the 
system  response  matrix  element).  In  the  3D  acquisition  case 
3  dimensions  are  associated  with  7,  while  4  dimensions  are 
associated  with  /  (as  can  readily  be  seen  for  one  2D  plane  in 
coicidence  with  another  2D  plane). 

The  EM  algorithm  generates  a  sequence  of  estimates  for 
each  iteration,  with  the  iterative  step  given  by: 

V‘+1)  =  o) 

Ar<*> 

r=  1 

The  iterative  step  for  the  ISRA  algorithm  can  be  written  as: 
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In  each  case  the  forward  projection  step  is  given  by 
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and  the  backprojection  step  is  given  by 
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In  Inverse  Monte  Carlo  sampling,  we  are  able  to  avoid 
explicit  calculation  (either  on-the-fly  or  precalculated  and 
stored)  of  the  system  response  matrix  elements  but  substituting 
a  statistical  estimate  of  the  desired  quantities.  A  Monte  Carlo 


sampled  estimate  of  the  forward  projection  is  straightforward, 
and  is  produced  by  generating  events  at  random  positions 
within  each  source  voxel,  with  the  number  of  events  generated 
per  voxel  proportional  to  the  activity  in  each  voxel  of  the 
source  estimate.  The  generated  events  are  then  propagated 
to  where  they  are  either  detected  or  escape  detection  (i.e. 
’’hit-or-miss’*).  Note  that  a  sophisticated  system  model  could 
here  follow  the  fate  of  each  photon  pair  in  considerable 
detail,  including  potential  attenuation  and  scatter  effects. 
The  Monte  Carlo  sampled  estimate  of  the  back-projection 
is  also  quite  straightforward:  a  fixed  number  of  events  are 
generated  randomly  within  each  source  voxel  in  turn,  with  the 
back-projection  amplitude  at  the  source  voxel  accumulated 
according  to  the  amplitude  of  the  line-of-response  bin  where 
each  successful  photon  pair  is  detected.  Note  that  the  system 
response  matrix  element  is  encoded  implicitly  in  the  probability 
that  a  Monte  Carlo  event  from  source  voxel  j  is  successfully 
detected  within  projection  i. 

Strictly  speaking,  the  above  description  of  the  EM  and 
ISRA  algorithms  assume  that  all  photons  produced  are  in 
fact  detected  at  some  line-of-response.  When  Monte  Carlo 
sampling,  both  the  forward-  and  back-projection  steps 
have  arbitrary  overall  normalization,  since  only  the  relative 
probabilities  are  important.  Therefore,  one  can  simply  set 
the  overall  number  of  forward-projected  events  in  the  4D 
space  equal  to  the  number  of  source  events  in  the  3D  space 
after  generating  any  desired  number  of  forward  projection 
sample  events.  Similarly,  one  can  generate  any  desired 
number  of  events  in  the  back-projection  step  and  subsequently 
preserve  normalization  by  setting  the  number  of  events  in  the 
back-projected  3D  space  equal  to  the  number  of  events  in  the 
4D  space  from  which  they  were  back-projected.  If  one  follows 
this  procedure,  the  iteration  steps  given  above  are  valid. 

III.  Methods 

We  have  been  evaluating  the  performance  of  Monte  Carlo 
estimated  forward-  and  back-projection  operations  with  explicit 
matrix  multiplication  using  pre-calculated  system  matrix 
elements.  These  matrix  elements  are  calculated  using  the  same 
Monte  Carlo  procedure  which  is  used  in  the  Monte  Carlo 
estimates,  and  at  the  moment  include  3D  interaction  points  for 
each  of  the  detected  gamma  photons  (i.e.  crystal  attenuation 
and  depth-of-interaction  effects  are  included)  but  with  only 
photocapture  and  no  detector  Compton  scatter  effects.  Events 
are  accumulated  in  all  4D  lines  of  response  for  a  given  voxel, 
but  the  results  are  stored  only  for  those  lines-of-response 
which  have  non-zero  values.  When  one  neglects  detector 
Compton  scatter  effects  (but  includes  inter-crystal  penetration 
and  depth-of-interaction  effects)  this  decreases  the  number  of 
stored  lines-of-response  dramatically. 

The  system  response  is  calculated  for  a  representative 
column  of  2  x  2  x  20  voxels,  each  1.5mm  x  1.5mm  x  1.5mm 
and  located  near  the  center  of  a  virtual  36  x  36  crystal  array 
of  3mm  x  3mm  x  15mm  crystals.  This  virtual  crystal  array  is 
4  times  the  size  of  our  actual  crystal  array,  allowing  a  single 
representative  voxel  at  each  depth  to  store  information  on  the 
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system  response  matrix  element  for  any  crystal  at  its  depth  in 
the  array  (e.g.  along  the  40-voxel-long  axis  between  the  two 
detector  planes).  Vertical  reflection  symmetry  is  used  to  reduce 
the  number  of  required  elements  by  two.  Translation  symmetry 
is  then  used,  with  each  matrix  element  checked  to  see  whether 
it  is  still  non-zero  after  translation  to  the  target  voxel  location. 
This  last  just  means  checking  whether  the  translated  photon 
pair  directions  (associated  with  a  given  matrix  element)  each 
still  hit  the  two  18x18  element  detector  arrays  after  translating 
the  source  voxel  to  a  given  location.  In  practice,  one  is  able 
to  calculate  what  range  of  translations  in  each  direction  will 
satisfy  this  condition,  obviating  the  need  for  a  explicit  check 
after  each  translation  and  thereby  speeding  the  algorithm 
execution. 

Having  implemented  the  above  procedures  and  done  some 
preliminary  qualitative  work  to  confirm  their  rough  validity, 
we  are  now  prepared  to  begin  quantitative  measurement  of  the 
relative  speed  and  performance  of  each  method.  This  will  be 
particularly  interesting  as  a  function  of  the  number  of  collected 
events,  and  as  a  function  of  the  accuracy  of  the  system  model 
incorporated  in  each  of  the  two  approaches.  To  carry  out  this 
study  we  plan  to  implement  procedures  analogous  to  those 
others  have  used  to  evaluate  3D  reconstruction  algorithms 
applied  to  larger  PET  systems,  using  reasonably  realistic 
software  phantoms  with  a  range  of  sizes,  orientations,  and 
contrasts.  In  particular,  we  are  implementing  a  series  of 
quantitative  estimators  of  hot-spot  detectability,  cold-spot 
detectabilty,  and  structural  accuracy  [10]. 

IV.  Discussion 

We  have  implemented  reconstruction  algorithms  for  a 
pair  of  18  x  18  crystal  arrays  used  in  a  fixed  parallel  planar 
acquisition  geometry.  Note  that  the  approach  using  explicit 
pre-calculation  of  the  system  matrix  elements  might  also 
be  used  for  larger  systems  if  one  preserves  the  translational 
symmetry.  In  practice,  this  means  limiting  accepted  events 
to  have  lines-of-response  with  the  top  crystal  no  more  than 
18  crystals  displaced  from  the  bottom  crystal  in  either  lateral 
dimension.  This  allows  for  ample  parallax  to  permit  accurate 
imaging  in  the  depth-of-field,  while  reducing  randoms  and 
scatter  acceptance  in  a  larger  field-of-view  device.  We  are 
currently  implementing  a  reconfigurable  array  of  18  x  18 
crystal  modules,  and  we  will  explore  application  of  these 
methods  to  this  larger  field-of-view  device. 

.  More  generally,  application  of  our  method  of  drastically 
compressing  the  system  response  matrix  depends  upon  the 
translational  invariance  (up  to  detector  edge  effects)  of  the 
system  response  matrix.  The  Monte  Carlo  sampling  approach 
is  not  so  constrained,  and  readily  incorporates  a  non-repeating 
and  largely  non-zero  system  response  matrix.  Such  a  matrix 
arises  when  including  such  effects  as  source-dependent 
attenuation  and  scatter.  For  this  reason,  we  are  particularly 
interested  in  evaluating  the  relative  capabilities  of  Monte  Carlo 
sampled  and  explicitly  calculated  versions  of  3D  iterative 
reconstruction  algorithms  within  this  simple,  constrained 
geometry. 


We  have  implemented  3D  iterative  reconstruction 
procedures  for  a  simple  fixed  dual-plane  acquisition  geometry, 
both  using  explicit  pre-calculation  of  system  response  matrix 
elements  and  using  Monte  Carlo  sampling  techniques.  The 
relative  speed  and  performance  of  these  techniques  will 
determine  their  range  of  application  as  a  function  of  event 
statistics  and  intrinsic  detector  system  resolution. 
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Abstract 

We  have  developed  a  three  dimensional  (3D)  reconstruction 
procedure  for  positron  emission  tomography  (PET)  based  on 
Inverse  Monte  Carlo  analysis  (IMC).  This  procedure  uses 
Monte  Carlo  techniques  to  simulate  the  process  of  photon 
emission  and  detection  from  a  hypothesized  source  distribution. 
A  Monte  Carlo  simulation  of  a  given  source/detector  system 
implicitly  contains  the  information  in  the  system  response 
matrix.  By  generating  Monte  Carlo  events  one  can  sample 
from  the  system  matrix  to  estimate  forward  or  backward 
projection  operations.  We  have  explored  this  approach  in  the 
context  of  the  image  space  reconstruction  algorithm  (ISRA)  [2] 
where  the  system  matrix  is  too  large  to  be  calculated  and  stored 
practically. 

This  technique  was  tested  on  a  series  of  software 
phantoms  and  also  from  data  taken  from  high  resolution 
PET  scanners  [1].  The  current  work  involves  quantitatively 
evaluating  the  performance  of  the  technique,  using  the 
methodology  suggested  by  Furie  et  al  [3,  4].  Reconstruction 
parameters,  like  the  number  of  Monte  Carlo  events  to  be 
generated  per  iterations,the  initial  source  hypothesis,  and  the 
stopping  criteria,  will  be  optimized  using  this  methodology. 
Performance  of  the  technique  will  be  evaluated  according  to 
figures  of  merit  associated  with  structural  accuracy,  hot  spot 
detectability  and  cold  spot  detectability. 

I.  INTRODUCTION 

The  recent  trend  in  operating  PET  detectors  in  three- 
dimensional  mode,  to  achieve  higher  sensitivity,  has 
heightened  interest  in  3D  image  reconstruction  algorithms. 
The  requirements  for  such  algorithms  is  that  they  have  to  be 
quantitatively  accurate,  and  have  fast  reconstruction  times. 
The  3DRP  algorithm  proposed  by  Kinahan  and  Rogers 
[5]  is  currently  commonly  implemented  on  commercial 
PET  scanners.  This  algorithm  is  fast,  robust  and  has  good 
image  contrast,  but  suffers  from  poor  noise  properties  and 
treatment  of  statistical  data.  Statistical  methods,  like  the 
maximum  likelihood  -  expectation  maximization  (ML-EM) 
[6,  7,  8],  penalized  weighted  least  squares  (PWLS)  [9]  and 
Bayesian  reconstruction  procedures [10],  give  improvements  in 
image  quality  at  increased  computational  cost.  These  set  of 
algorithms  are  iterative,  and  require  explicit  calculation  of  the 
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system  response  matrix.  This  matrix  gives  the  probability  that 
a  photon  emitted  from  a  certain  source  element  will  be  detected 
in  a  particular  detector  line  of  response.  This  probability 
matrix  is  usually  much  too  large  to  store  in  memory  and  has 
to  be  either  stored  on  disk  or  calculated  on  the  fly.  Several 
groups  have  used  symmetries  and  sparseness  of  this  matrix 
to  be  able  to  store  it  in  memory  [10].  The  ordered  subsets 
(OS-EM)  algorithm  [11]  uses  subsets  of  the  data  to  speed  up 
the  iterations  but  the  reconstruction  times  for  3D  mode  are  still 
quite  large.  Approximate  rebinning  algorithms  like  the  FORE 
algorithm  rebin  data  into  set  of  2D  slices  [12];  the  rebinned 
slices  can  be  reconstructed  using  filtered  backprojection  or 
ML-EM,  OS-EM  or  WLS  methods.  These  algorithms  are 
computationally  fast  but  the  accuracy  of  the  reconstructed 
images  is  limited  by  the  approximations  implicit  in  the  line 
integral  model  and  the  non -Poisson  statistical  properties  of  the 
resultant  slice  data  [13]. 

Our  interest  in  image  reconstruction  arises  from 
development  of  a  novel  detector  with  dual  plane  acquisition 
geometry  [14].  For  this  detector,  the  incomplete  set  of 
projections  excludes  the  possibility  of  using  Fourier  filtered 
backprojection  techniques  since  the  geometry  does  not 
satisfy  Orlov’s  criterion  [15].  Iterative  algorithms  which  are 
not  restricted  by  device  geometry  are  a  possibility  but  the 
computational  burden  associated  with  explicitly  calculating 
the  .  system  response  function  makes  them  difficult  if  not 
impractical. 

Monte  Carlo  methods  for  solving  the  inverse  problem  in 
ECT  has  been  previously  reported  by  other  groups  [16].  This 
method  generally  uses  Monte  Carlo  methods  to  calculate  the 
system  response  matrix,  taking  into  account  attenuation  and 
scatter  of  the  gamma  rays.  The  ML-EM  algorithm  is  then  used 
to  solve  for  the  source  distribution.  Monte-Carlo  methods  have 
also  been  used  in  PET  to  study  the  effects  of  attenuation  and 
scatter  as  they  provide  a  natural  framework  for  incorporation  of 
detector  response  functions,  and  the  effects  of  attenuation  and 
other  sources  of  systematic  errors. 

Our  approach  may  be  distinguished  from  these  approaches, 
in  that,  it  uses  Monte  Carlo  generation  of  events  to  model  the 
emission  and  detection  of  photons  from  a  hypothesized  source 
distribution  within  each  iteration,  and  thus  avoids  the  explicit 
calculation  of  the  system  response  matrix. 

II.  Theory 

Our  first  application  of  Inverse  Monte  Carlo  procedures  has 
been  in  the  context  of  the  3D  ISRA  algorithm.  This  algorithm 
introduced  the  concept  of  storing  data  as  backprojection 
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array,  thereby  compressing  the  large  and  sometimes  sparsely 
populated  projection  matrix  into  a  more  compact  data  image 
matrix.  In  this  algorithm,  the  update  to  the  source  hypothesis 
is  done  by  comparing,  the  backprojection  of  the  data  and  the 
backprojection  of  projections  from  the  source  hypothesis. 
Titterington  [17]  showed  that  this  algorithm  converges  to  the 
least  squares  estimate,  as  long  as  the  source  elements  are 
strictly  positive.  Using  his  notation;  for  J  source  pixels,  the 
jth  element  has  emission  density  A  j.  The  measured  data  nf 
is  the  number  of  coincidences  counted  in  the  ith  of  7  pairs  of 
coincident  detector  elements.  The  conditional  probability  that 
an  event  emitted  from  pixel  j  is  assigned  to  projection  i  is  given 
by  aij.  The  iterative  step  for  the  ISRA  algorithm  can  be  written 
as 
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The  IMC  technique  replaces  the  explicit  forward  projection 
calculation  in  ISRA,  by  Monte  Carlo  generation  of  events 
according  to  the  hypothesized  source  distribution.  Photon 
pairs  are  generated  at  random  positions  within  the  source 
voxels,  with  the  number  of  photon  pairs  generated  per  voxel 
proportional  to  the  hypothesized  source  intensity  at  that 
voxel. These  photon  pairs  are  given  a  random  direction  and 
traced  to  determine  where  and  whether  they  hit  the  detector. 
An  elaborate  source  and  system  response  model  could  contain 
attenuation  and  scatter  in  the  source,  as  well  as  details  of 
the  detector  response  properties  like  resolution,  scatter  and 
attenuation  in  the  detector  elements.  In  practice,  we  have 
explored  a  very  simple  system  model  containing  only  the  basic 
detector  geometry. 

The  events  that  hit  the  detector  are  backprojected.  In  formal 
terms  backprojection  is  the  convolution  of  the  system  response 
matrix  with  the  projection  matrix  as  shown  in  equation  3.  The 
system  response  matrix  element  a,j  can  be  approximated  by  the 
length  of  the  intersection  of  the  ith  line  with  the  jth  pixel  [18]. 
Hence  the  backprojection  operation  can  be  approximated  as,  the 
accumulation  according  to  the  chord  length  by  which  each  cubic 
voxel  along  a  line  of  response  is  intercepted.  [19]. 

Like  in  ISRA,  the  ratio  the  backprojection  of  the  data  and 
the  backprojection  of  the  simulated  data,  for  a  voxel  element,  is 
used  to  update  the  image  voxel  element. 

The  backprojection  of  the  data  into  image  space  is  used 
as  the  initial  source  estimate  during  the  first  iteration.  In  the 
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following  iterations,  the  number  of  events  generated  from 
a  voxel  is  updated  (multiplied)  by,  the  ratio  between  the 
backprojection  of  the  data  to  the  backprojection  of  the  Monte 
Carlo  generated  “pseudo-data”at  that  voxel.  The  pseudo-data 
events  are  similarly  projected  to  find  where  they  hit  the  detector 
and  then  backprojected.If  the  simulated  backprojection  at  the 
source  voxel  is  less  than  the  measured  backprojection  at  that 
voxel,  then  more  pseudo-data  events  are  generated.  Conversely, 
if  the  simulated  backprojection  is  more  than  the  measured 
backprojection  then,  the  pseudo-data  events  are  removed. 

III.  Preliminary  Results 

A  simulated  Hoffman  brain  phantom  [20]  with  1  mm 
spatial  resolution;  neglecting  attenuation  and  scatter,  was 
reconstructed  using  the  IMC  algorithm.  The  detector  geometry 
used  for  this  purpose  was  that  of  the  HEAD  PENN-PET 
detector.  The  HEAD  PENN-PET  scanner  is  a  volume-imaging 
PET  scanner  without  inter-plane  septa,  with  high  spatial 
resolution  and  sensitivity  [21].  75  million  events  were 

simulated  from  the  phantom  and  the  data  was  stored  as  a 
backprojection  array.  The  images  shown  in  Figure  1  were 
reconstructed  simulating  350  million  events,  using  4  mm3 
voxel  size,  with  the  image  array  size  of  60  X  60  X  48.  The 
algorithm  was  implemented  on  a  SGI/Cray  Grigin2  Series 
parallel  computer  using  8  processors.  The  reconstruction  time 
required  for  this  implementation  was  about  5  minutes  for  25 
iterations  . 


Figure  1:  Simulated  Hoffman  Brain  Phantom  using  4  mm3  voxel  size 
and  60  X  60  X  48  image  array  size  ;  using  25  iterations  and  simulation 
350  million  events 


The  algorithm  was  also  applied  to  Hoffman  brain  data 
acquired  by  the  HEAD  PENN-PET  data.  The  total  number 
of  attenuation  corrected  events  in  this  data  were  about  800 
million.  Using  4  mm  size  voxels  the  image  (60  X60  48) 
was  reconstructed  in  25  iterations  within  15  minutes  while 
simulating  3  billion  events  [1].  The  images  from  this 
reconstruction  are  shown  in  Figure  2.  In  this  implementation 
the  simulated  events  are  traced  to  find  where  they  hit  the 


detector  and  accumulated.  Backprojection  is  done  once  per 
iteration  after  all  the  simulated  events  are  generated.  Using 
symmetries  in  the  image  array  the  backprojection  is  sped  by  an 
order  of  magnitude.  Hence  reconstruction  time  does  not  scale 
up  linearly  with  increased  number  of  simulated  events. 


Figure  2:  Hoffman  Brain  Phantom  using  4  mm3  voxel  size  and 
60  X60  X48  array  image  size;  using  25  iterations  and  simulating  3 
billion  events.  Data  was  acquired  by  HEAD  PENN-PET  detector 


IV.  Discussion 

The  images  shown  above  demonstrate  the  feasibility  of 
using  Inverse  Monte  Carlo  for  image  reconstruction.  The  next 
task,  which  is  currently  underway,  is  to  quantitatively  evaluate 
the  performance  of  this  algorithm.  To  do  so  we  are  using 
the  methodology  proposed  by  the  Medical  Image  Processing 
Group  (MIPG)  at  the  University  of  Pennsylvania  [3,  4].  They 
have  developed  a  software  package  called  EVAL3DPET  [22] 
which  are  a  set  of  programs  designed  to  statistically  evaluate 
3D  PET  reconstruction  algorithms.  This  software  is  capable  of 
generating  phantom  and  projection  data  based  on  a  realistic  3D 
PET  scanner  model.  The  phantoms  are  random  samples  from  a 
statistically  described  ensemble  of  3D  images  with  69  ellipsoid 
features  with  varying  activity  and  sizes.  The  projection  data 
takes  into  account  detector  field-of-view  blurring  and  a  realistic 
3D  PET  noise  model. 

The  images  reconstructed  from  this  phantom  can  be 
evaluated  for  structural  accuracy,  hot  spot  detectabilty  and  cold 
spot  detectabilty.  The  evaluation  program  also  calculates  a 
training  figure-of-meriuhat  can  be  used  for  optimizing  the  free 
parameters  of  the  reconstruction  techniques.  The  parameters  of 
the  IMC  algorithm  that  we  would  like  to  optimize  are  the  total 
number  of  events  to  simulate,  the  number  of  events  to  simulate 
per  iteration,  the  number  of  iterations  or  the  stopping  criteria 
and  also,  the  choice  of  our  first  estimate. 

Increasing  the  number  of  events  generated  reduces 
statistical  fluctuations  of  the  simulated  data,  and  yields  better 
images.  Generating  a  large  number  of  events,  however, 
increases  the  reconstruction  time,  without  appreciably 


improving  image  quality.  We  have  worked  with  a  stopping 
criteria  that  the  number  of  events  generated  is  about  10  times 
the  number  of  events  detected.  We  have  observed  empirically 
that  the  image  quality  does  not  improve  significantly  for  more 
events  generated. 

There  are  various  choices  for  the  first  estimate  of  the  source 
distribution  ,  such  as  a  uniform  field  or  an  image  formed  from 
some  earlier  iteration  or  other  algorithm  .The  backprojection  of 
the  data  can  also  be  used  as  the  first  estimate,  and  it  has  been 
observed  that  this  is  a  good  starting  point.  Using  the  training 
figures-of-merit  we  will  form  the  basis  of  selecting  the  best 
starting  point  to  give  fast,  accurate  images  of  3D  PET  data. 

V.  Conclusion 

We  have  explored  the  possibility  using  Inverse  Monte 
Carlo  methods  for  image  reconstruction  and  applied  it  to 
digital  phantoms,  as  well  as,  data  from  3D  PET  detectors. 
The  algorithm  is  shown  to  give  qualitatively  accurate  images 
of  complex  phantoms.  This  work  will  present  quantitative 
evaluation  of  the  algorithm  using  figures  of  merit  based 
on  structural  accuracy,  hot  spot  detectabilty  and  cold  spot 
detectabilty. 
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CONCLUSIONS 


Our  principal  accomplishment  has  been  the  successful  construction  and  imaging 
characterization  of  a  12cm  x  12cm  field-of-view  PEM  instrument.  Performance  measurements  and 
other  outcomes  for  this  preliminary  device  are  detailed  below  by  category. 

1)  Spatial  Resolution  and  Uniformity: 

We  have  obtained  2mm  FWHM  (Full  Width  at  Half  Maximum)  spatial  resolution  in 
imaging  a  point-like  source.  Quantitatively  uniform  sensitivity  over  a  large  spatial  range  was  also. 
Spatial  calibration  of  the  detector  was  straightforward  when  imaging  a  point  source,  and  3-4mm 
FWHM  resolution  was  obtainable  with  no  calibration  constants  at  all.  Boundary  effects  between 
adjacent  fiber  ribbons  were  found  to  be  minimal,  and  determination  of  interaction  coordinates  was 
robust  and  unambiguous  because  of  the  substantial  fiber  light  yield  (averaging  8 
photoelectrons/511  keV  interaction)  and  because  of  the  intrinsic  linearity  of  interaction  position 
measurement  through  local  light  collection  by  fibers.  The  coordinates  obtained  were  also  accurate 
to  the  very  edge  of  each  prototype  detector,  as  opposed  to  the  Anger  coordinates  which  were  quite 
distorted  near  each  edge. 

2)  Light  Yield: 

We  have  obtained  an  average  of  8  photoelectrons  at  the  fiber  readout  for  1mm  fibers  and  roughly 
200  photoelectrons  at  the  Anger  PMT  readout,  using  CsI(Na)  crystals.  The  Anger  light  yield  was 
low  partially  because  of  incomplete  photocathode  coverage  by  four  2”  round  PMTs,  and  we  have 
purchased  and  obtained  60mm  square  PMTs  for  successor  modules.  Another  cause  is  the 
relatively  low  PMT  quantum  efficiency  for  green  wave-shifted  light. 

3)  Energy  Resolution: 

We  have  obtained  20%  energy  resolution  FWHM,  after  correction  for  nonuniformity  in  the 
Anger  PMT  light  collection  across  the  module  (this  calibration  is  straightforward,  given  the  fiber- 
derived  event  coordinates  and  the  511  keV  photopeak).  Of  this,  about  15%  FWHM  is  intrinsic 
nonuniformity  and  about  15%  is  due  to  our  limited  photostatistics  (these  effects  add  in  quadrature). 
Therefore,  we  anticipate  a  future  energy  resolution  of  about  15%  FWHM,  assuming  400 
photoelectrons/5 1 1  keV  collection.  This  should  easily  be  adequate  for  clinical  scatter  rejection. 

4)  Timing  Resolution: 

We  have  obtained  50ns  FWHM  timing  resolution  for  CsI(Na)  modules  by  using  a  second  low 
threshold  on  the  Anger  sums  in  addition  to  the  threshold  used  for  energy ’trigger  discrimination. 
This  is  in  agreement  with  the  decay  time  of  the  CsI(Na)  and  the  results  from  computer  simulation 
of  the  best  timing  achievable  as  a  function  of  integration  time  and  threshold  level. 

5)  Optoelectronics  and  Data  Acquisition  Electronics: 

We  have  developed  a  high-performance,  low-cost,  and  compact  data  acquisition  system  suitable 
for  high-rate  PET  data  acquisition  with  a  portable  instrument. 

6)  Image  Reconstruction  and  Analysis: 

Reconstruction  algorithms,  including  background  subtraction  and  attenuation/scatter 
corrections,  is  being  carried  out  with  a  novel  iterative  reconstruction  algorithm  which  we  have 
developed.  Iterative  reconstruction  permits  3d  limited-angle  tomography,  such  as  is  our  case  with 
two  modules  which  do  not  completely  surround  the  region  of  interest  (i.e.  a  breast  or  axillary 
tails). 
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